Question

求用C语言实现FFT变换的程序（见下面）

2.用FFT程序分析正 玄信号，分别在以下情况进行分析，并讨论所得的结果
a 信号频率F＝50Hz，采样长N＝32,采样时间T＝0.000625s
b 信号频率F＝50Hz，采样长N＝32,采样时间T＝0.005s
c 信号频率F＝50Hz，采样长N＝32,采样时间T＝0.0046875s
d 信号频率F＝50Hz，采样长N＝32,采样时间T＝0.004s
e 信号频率F＝50Hz，采样长N＝64,采样时间T＝0.000625s
f 信号频率F＝250Hz，采样长N＝32,采样时间T＝0.005s
g 将C信号后补32个0,做64点FFT
三 实验报告
1.记录下实验内容中各种情况下的X(k)值，作出频谱图，并深入讨论结果，说明参数的变化对信号频谱产生哪些影响。频谱只作模特性，模的最大值＝1,全部规一化。
2.打印出用C语言编写的FFT源程序，并且在每一小段处加上详细的注释说明。
3.用C语言编写FFT编写程序时，要求采用人机界面形式：
N、T、F变量均由键盘输入；补零或不补零要求设置一开关。

Answer 1

你好，这是我的回答，希望可以帮到你。

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教你。

Answer 2

这是一个傅里叶变化的子函数，你可以自己做主函数传递你这里的参数验证
// 入口参数：
// l: l = 0, 傅立叶变换; l = 1, 逆傅立叶变换
// il: il = 0,不计算傅立叶变换或逆变换模和幅角；il = 1,计算模和幅角
/歼码旦/ n: 输入的点数，为偶数，一般为32，64，128，...,1024等
// k: 满足n=2^k(k>0),实质上k是n个采样数据可以分解为偶次幂和奇次幂的次数
// pr[]: l=0时，存放N点采样数据的实部
// l=1时, 存放傅立叶变换的N个实部
// pi[]: l=0时，存放N点采样数据的虚部
// l=1时, 存放傅立叶变换的N个虚部
//
// 出口参数：
// fr[]: l=0, 返回傅立叶变换的实部
// l=1, 返回逆傅立叶变换的实部
// fi[]: l=0, 返回傅立叶变换的虚部
// l=1, 返回逆傅立叶变换的虚部
// pr[]: il = 1,l = 0 时，返回傅立叶变换的模
// il = 1,l = 1 时，返回逆傅立叶变换的模
// pi[]: il = 1,l = 0 时，返回傅立叶变换的辐角
// il = 1,l = 1 时，返回逆傅立氏扰叶变换的辐角

void kbfft(double *pr,double *pi,int n,int k,double *fr,double *fi,int l,int il)
{
int it,m,is,i,j,nv,l0;
double p,q,s,vr,vi,poddr,poddi;

//排序
for (it=0; it<=n-1; it++)
{ m=it; is=0;
for (i=0; i<=k-1; i++)
{ j=m/2; is=2*is+(m-2*j); m=j;
fr[it]=pr[is]; fi[it]=pi[is];
}
}

//模差蝶形运算
pr[0]=1.0; pi[0]=0.0;
p=6.283185306/(1.0*n);
pr[1]=cos(p); pi[1]=-sin(p);
if (l!=0) pi[1]=-pi[1];
for (i=2; i<=n-1; i++)
{ p=pr[i-1]*pr[1]; q=pi[i-1]*pi[1];
s=(pr[i-1]+pi[i-1])*(pr[1]+pi[1]);
pr[i]=p-q; pi[i]=s-p-q;
}
for (it=0; it<=n-2; it=it+2)
{ vr=fr[it]; vi=fi[it];
fr[it]=vr+fr[it+1]; fi[it]=vi+fi[it+1];
fr[it+1]=vr-fr[it+1]; fi[it+1]=vi-fi[it+1];
}
m=n/2; nv=2;
for (l0=k-2; l0>=0; l0--)
{ m=m/2; nv=2*nv;
for (it=0; it<=(m-1)*nv; it=it+nv)
for (j=0; j<=(nv/2)-1; j++)
{ p=pr[m*j]*fr[it+j+nv/2];
q=pi[m*j]*fi[it+j+nv/2];
s=pr[m*j]+pi[m*j];
s=s*(fr[it+j+nv/2]+fi[it+j+nv/2]);
poddr=p-q; poddi=s-p-q;
fr[it+j+nv/2]=fr[it+j]-poddr;
fi[it+j+nv/2]=fi[it+j]-poddi;
fr[it+j]=fr[it+j]+poddr;
fi[it+j]=fi[it+j]+poddi;
}
}
if (l!=0)
for (i=0; i<=n-1; i++)
{ fr[i]=fr[i]/(1.0*n);
fi[i]=fi[i]/(1.0*n);
}
if (il!=0)
for (i=0; i<=n-1; i++)
{ pr[i]=sqrt(fr[i]*fr[i]+fi[i]*fi[i]);
pr[i]=(pr[i]/(n/2)); //各次谐波幅值，其中pr[1]为基波幅值
if (fabs(fr[i])<0.000001*fabs(fi[i]))//fabs()是取绝对值函数,浮点型的0 在内存中并不是严格等于0,可以认为当一个浮点数

离原点足够近时,也就是f>0.00001 && f<-0.00001,认为f是0
{ if ((fi[i]*fr[i])>0) pi[i]=90.0;
else pi[i]=-90.0;
}
else
pi[i]=atan(fi[i]/fr[i])*360.0/6.283185306;
}
return;
}