扫二维码与项目经理沟通
我们在微信上24小时期待你的声音
解答本文疑问/技术咨询/运营咨询/技术建议/互联网交流
你好,这是我的回答,希望可以帮到你。
创新互联长期为上千家客户提供的网站建设服务,团队从业经验10年,关注不同地域、不同群体,并针对不同对象提供差异化的产品和服务;打造开放共赢平台,与合作伙伴共同营造健康的互联网生态环境。为夏邑企业提供专业的成都做网站、网站建设、外贸营销网站建设,夏邑网站改版等技术服务。拥有十载丰富建站经验和众多成功案例,为您定制开发。
1)结果讨论
一,如果对信号进行同样点数N的FFT变换,采样频率fs越高,则可以分析越高频的信号;与此同时,采样频率越低,对于低频信号的频谱分辨率则越好。
二,假设采样点不在正弦信号的波峰、波谷、以及0电压处,频谱则会产生泄露(leakage)。
三,对于同样的采样率fs,提高FFT的点数N,则可提高频谱的分辨率。
四,如果采样频率fs小于2倍信号频率2*fs(奈圭斯特定理),则频谱分析结果会出错。
五,对于(二)中泄露现象,可以通过在信号后面补零点解决。
2)程序及注解如下
%清除命令窗口及变量
clc;
clear all;
%输入f、N、T、是否补零(补几个零)
f=input('Input frequency of the signal: f\n');
N=input('Input number of pointsl: N\n');
T=input('Input sampling time: T\n');
flag=input('Add zero too sampling signal or not? yes=1 no=0\n');
if(flag)
ZeroNum=input('Input nmber of zeros\n');
else
ZeroNum=0;
end
%生成信号,signal是原信号。signal为采样信号。
fs=1/T;
t=0:0.00001:T*(N+ZeroNum-1);
signal=sin(2*pi*f*t);
t2=0:T:T*(N+ZeroNum-1);
signal2=sin(2*pi*f*t2);
if (flag)
signal2=[signal2 zeros(1, ZeroNum)];
end
%画出原信号及采样信号。
figure;
subplot(2,1,1);
plot(t,signal);
xlabel('Time(s)');
ylabel('Amplitude(volt)');
title('Singnal');
hold on;
subplot(2,1,1);
stem(t2,signal2,'r');
axis([0 T*(N+ZeroNum) -1 1]);
%作FFT变换,计算其幅值,归一化处理,并画出频谱。
Y = fft(signal2,N);
Pyy = Y.* conj(Y) ;
Pyy=(Pyy/sum(Pyy))*2;
f=0:fs/(N-1):fs/2;4
subplot(2,1,2);
bar(f,Pyy(1:N/2));
xlabel('Frequency(Hz)');
ylabel('Amplitude');
title('Frequency compnents of signal');
axis([0 fs/2 0 ceil(max(Pyy))])
grid on;
祝你好运!
我可以帮助你,你先设置我最佳答案后,我百度Hii教你。
1、二维FFT相当于对行和列分别进行一维FFT运算。具体的实现办法如下:
先对各行逐一进行一维FFT,然后再对变换后的新矩阵的各列逐一进行一维FFT。相应的伪代码如下所示:
for (int i=0; iM; i++)
FFT_1D(ROW[i],N);
for (int j=0; jN; j++)
FFT_1D(COL[j],M);
其中,ROW[i]表示矩阵的第i行。注意这只是一个简单的记法,并不能完全照抄。还需要通过一些语句来生成各行的数据。同理,COL[i]是对矩阵的第i列的一种简单表示方法。
所以,关键是一维FFT算法的实现。
2、例程:
#include stdio.h
#include math.h
#include stdlib.h
#define N 1000
/*定义复数类型*/
typedef struct{
double real;
double img;
}complex;
complex x[N], *W; /*输入序列,变换核*/
int size_x=0; /*输入序列的大小,在本程序中仅限2的次幂*/
double PI; /*圆周率*/
void fft(); /*快速傅里叶变换*/
void initW(); /*初始化变换核*/
void change(); /*变址*/
void add(complex ,complex ,complex *); /*复数加法*/
void mul(complex ,complex ,complex *); /*复数乘法*/
void sub(complex ,complex ,complex *); /*复数减法*/
void output();
int main(){
int i; /*输出结果*/
system("cls");
PI=atan(1)*4;
printf("Please input the size of x:\n");
scanf("%d",size_x);
printf("Please input the data in x[N]:\n");
for(i=0;isize_x;i++)
scanf("%lf%lf",x[i].real,x[i].img);
initW();
fft();
output();
return 0;
}
/*快速傅里叶变换*/
void fft(){
int i=0,j=0,k=0,l=0;
complex up,down,product;
change();
for(i=0;i log(size_x)/log(2) ;i++){ /*一级蝶形运算*/
l=1i;
for(j=0;jsize_x;j+= 2*l ){ /*一组蝶形运算*/
for(k=0;kl;k++){ /*一个蝶形运算*/
mul(x[j+k+l],W[size_x*k/2/l],product);
add(x[j+k],product,up);
sub(x[j+k],product,down);
x[j+k]=up;
x[j+k+l]=down;
}
}
}
}
/*初始化变换核*/
void initW(){
int i;
W=(complex *)malloc(sizeof(complex) * size_x);
for(i=0;isize_x;i++){
W[i].real=cos(2*PI/size_x*i);
W[i].img=-1*sin(2*PI/size_x*i);
}
}
/*变址计算,将x(n)码位倒置*/
void change(){
complex temp;
unsigned short i=0,j=0,k=0;
double t;
for(i=0;isize_x;i++){
k=i;j=0;
t=(log(size_x)/log(2));
while( (t--)0 ){
j=j1;
j|=(k 1);
k=k1;
}
if(ji){
temp=x[i];
x[i]=x[j];
x[j]=temp;
}
}
}
/*输出傅里叶变换的结果*/
void output(){
int i;
printf("The result are as follows\n");
for(i=0;isize_x;i++){
printf("%.4f",x[i].real);
if(x[i].img=0.0001)printf("+%.4fj\n",x[i].img);
else if(fabs(x[i].img)0.0001)printf("\n");
else printf("%.4fj\n",x[i].img);
}
}
void add(complex a,complex b,complex *c){
c-real=a.real+b.real;
c-img=a.img+b.img;
}
void mul(complex a,complex b,complex *c){
c-real=a.real*b.real - a.img*b.img;
c-img=a.real*b.img + a.img*b.real;
}
void sub(complex a,complex b,complex *c){
c-real=a.real-b.real;
c-img=a.img-b.img;
}
ELEM这个数组是为了确定二维数组的每个值的确定位置,因为你申请的in和out都是1维数组(不考虑实部虚部那个0和1)所以要确定在你开辟的这个1维数组中所对应的二维数组是哪个数就用到了ELEM,举个例子,比如你算二维数组[r,c]这个位置的值,其实是在一维数组的第r*N+c个数,所以就得到了in[r*N+c][0],也就是在做fft之前的fftw_complex表示,楼主可能混淆fftw_complex和文中ELEM数组了,ELEM并不是fftw_complex数组的表示,只是为了找到确定数组位置的一个中间数组而已,希望回答能让你满意
#include stdio.h
#include math.h
#include stdlib.h
#define N 1000
/*定义复数类型*/
typedef struct{
double real;
double img;
}complex;
complex x[N], *W; /*输入序列,变换核*/
int size_x=0; /*输入序列的大小,在本程序中仅限2的次幂*/
double PI; /*圆周率*/
void fft(); /*快速傅里叶变换*/
void initW(); /*初始化变换核*/
void change(); /*变址*/
void add(complex ,complex ,complex *); /*复数加法*/
void mul(complex ,complex ,complex *); /*复数乘法*/
void sub(complex ,complex ,complex *); /*复数减法*/
void output();
int main(){
int i; /*输出结果*/
system("cls");
PI=atan(1)*4;
printf("Please input the size of x:\n");
scanf("%d",size_x);
printf("Please input the data in x[N]:\n");
for(i=0;isize_x;i++)
scanf("%lf%lf",x[i].real,x[i].img);
initW();
fft();
output();
return 0;
}
/*快速傅里叶变换*/
void fft(){
int i=0,j=0,k=0,l=0;
complex up,down,product;
change();
for(i=0;i log(size_x)/log(2) ;i++){ /*一级蝶形运算*/
l=1i;
for(j=0;jsize_x;j+= 2*l ){ /*一组蝶形运算*/
for(k=0;kl;k++){ /*一个蝶形运算*/
mul(x[j+k+l],W[size_x*k/2/l],product);
add(x[j+k],product,up);
sub(x[j+k],product,down);
x[j+k]=up;
x[j+k+l]=down;
}
}
}
}
/*初始化变换核*/
void initW(){
int i;
W=(complex *)malloc(sizeof(complex) * size_x);
for(i=0;isize_x;i++){
W[i].real=cos(2*PI/size_x*i);
W[i].img=-1*sin(2*PI/size_x*i);
}
}
/*变址计算,将x(n)码位倒置*/
void change(){
complex temp;
unsigned short i=0,j=0,k=0;
double t;
for(i=0;isize_x;i++){
k=i;j=0;
t=(log(size_x)/log(2));
while( (t--)0 ){
j=j1;
j|=(k 1);
k=k1;
}
if(ji){
temp=x[i];
x[i]=x[j];
x[j]=temp;
}
}
}
/*输出傅里叶变换的结果*/
void output(){
int i;
printf("The result are as follows\n");
for(i=0;isize_x;i++){
printf("%.4f",x[i].real);
if(x[i].img=0.0001)printf("+%.4fj\n",x[i].img);
else if(fabs(x[i].img)0.0001)printf("\n");
else printf("%.4fj\n",x[i].img);
}
}
void add(complex a,complex b,complex *c){
c-real=a.real+b.real;
c-img=a.img+b.img;
}
void mul(complex a,complex b,complex *c){
c-real=a.real*b.real - a.img*b.img;
c-img=a.real*b.img + a.img*b.real;
}
void sub(complex a,complex b,complex *c){
c-real=a.real-b.real;
c-img=a.img-b.img;
}
我们在微信上24小时期待你的声音
解答本文疑问/技术咨询/运营咨询/技术建议/互联网交流