求用C语言实现FFT变换的程序(见下面)
Posted
tags:
篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了求用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变量均由键盘输入;补零或不补零要求设置一开关。
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教你。 参考技术A 这是一个傅里叶变化的子函数,你可以自己做主函数传递你这里的参数验证
// 入口参数:
// 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;
参考资料:http://hi.baidu.com/xfcylyf/blog/item/bb4aaef3b3d307cd0a46e001.html
以上是关于求用C语言实现FFT变换的程序(见下面)的主要内容,如果未能解决你的问题,请参考以下文章
求用C语言编写一个简单的学生信息管理程序【 C++】网上之前有C 的好像