#include <stdio.h>
#define N 8
void kkfft(double pr[], double pi[], int n, int k, double fr[], double fi[], int l, int il)
void main()
{
double xr[N],xi[N],Yr[N],Yi[N],l=0,il=0
int i,j,n=N,k=3
for(i=0i<Ni++)
{
xr[i]=i
xi[i]=0
}
printf("------FFT------\n"旁此)
l=0
kkfft(xr,xi,n,k,Yr,Yi,l,il)
for(i=0i<Ni++)
{
printf("%-11lf + j* %-11lf\n",Yr[i],Yi[i])
}
printf("-----DFFT-------\n")
l=1
kkfft(Yr,Yi,n,k,xr,xi,l,il)
for(i=0i<Ni++)
{
printf("%-11lf + j* %-11lf\谨伍n"祥启或,xr[i],xi[i])
}
getch()
}
void kkfft(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=0it<=n-1it++)
{
m = it
is = 0
for(i=0i<=k-1i++)
{
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=2i<=n-1i++)
{
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=0it<=n-2it=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-2l0>=0l0--)
{
m = m/2
nv = 2*nv
for(it=0it<=(m-1)*nvit=it+nv)
for (j=0j<=(nv/2)-1j++)
{
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=0i<=n-1i++)
{
fr[i] = fr[i]/(1.0*n)
fi[i] = fi[i]/(1.0*n)
}
}
/*是否计算模和相角*/
if(il!=0)
{
for(i=0i<=n-1i++)
{
pr[i] = sqrt(fr[i]*fr[i]+fi[i]*fi[i])
if(fabs(fr[i])<0.000001*fabs(fi[i]))
{
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
}
void fft(){
int nn,n1,n2,i,j,k,l,m,s,l1
float ar[1024],ai[1024] // 实部 虚部
float a[2050]
float t1,t2,x,y
float w1,w2,u1,u2,z
float fsin[10]={0.000000,1.000000,0.707107,0.3826834,0.1950903,0.09801713,0.04906767,0.02454123,0.01227154,0.00613588,}// 优化
float fcos[10]={-1.000000,0.000000,0.7071068,0.9238796,0.9807853,0.99518472,0.99879545,0.9996988,0.9999247,0.9999812,}
nn=1024
s=10
n1=nn/2 n2=nn-1
j=1
for(i=1i<=nni++)
{
a[2*i]=ar[i-1]
a[2*i+1]=ai[i-1]
}
for(l=1l<n2l++)
{
if(l<j)
{
t1=a[2*j]
t2=a[2*j+1]
a[2*j]=a[2*l]
a[2*j+1]=a[2*l+1]
a[2*l]=t1
a[2*l+1]=t2
}
斗友 k=n1
则物 while (k<j)
{
j=j-k
k=k/2
}
j=j+k
}
for(i=1i<=si++)
{
u1=1
u2=0
m=(1<<i)
k=m>>1
w1=fcos[i-1]
w2=-fsin[i-1]
for(j=1j<=kj++)
{
for(l=jl<nnl=l+m)
{
孙销液 l1=l+k
t1=a[2*l1]*u1-a[2*l1+1]*u2
t2=a[2*l1]*u2+a[2*l1+1]*u1
a[2*l1]=a[2*l]-t1
a[2*l1+1]=a[2*l+1]-t2
a[2*l]=a[2*l]+t1
a[2*l+1]=a[2*l+1]+t2
}
z=u1*w1-u2*w2
u2=u1*w2+u2*w1
u1=z
}
}
for(i=1i<=nn/2i++)
{
ar[i]=a[2*i+2]/nn
ai[i]=-a[2*i+3]/nn
a[i]=4*sqrt(ar[i]*ar[i]+ai[i]*ai[i]) // 幅值
}
}
先给你个利用matlab中傅里叶变换进行函数频谱分析的程序。clf
fs=100N=128 %采样频率和数据点数空罩
n=0:N-1t=n/fs %时间序列
x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t)%信号
y=fft(x,N) %对信号进行快速Fourier变换
mag=abs(y)%求得Fourier变换后的振幅
f=n*fs/N %频率序列
subplot(2,2,1),plot(f,mag) %绘出随频率斗困闹变化的振幅
xlabel('频率/尺衫Hz')
ylabel('振幅')title('N=128')grid on
subplot(2,2,2),plot(f(1:N/2),mag(1:N/2))%绘出Nyquist频率之前随频率变化的振幅
xlabel('频率/Hz')
ylabel('振幅')title('N=128')grid on
%对信号采样数据为1024点的处理
fs=100N=1024n=0:N-1t=n/fs
x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t)%信号
y=fft(x,N) %对信号进行快速Fourier变换
mag=abs(y) %求取Fourier变换的振幅
f=n*fs/N
subplot(2,2,3),plot(f,mag)%绘出随频率变化的振幅
xlabel('频率/Hz')
ylabel('振幅')title('N=1024')grid on
subplot(2,2,4)
plot(f(1:N/2),mag(1:N/2))%绘出Nyquist频率之前随频率变化的振幅
xlabel('频率/Hz')
ylabel('振幅')title('N=1024')grid on
系统实现的意义和必要性：通过变换可以将原本的时域信号函数转化为频域，可以直观的观察到所采集信号函数的频域特征，更有利于进行信号分析
系统功能：通过傅里叶变换将信号函数的时域图形转化成频域图形，即将信号函数原本幅值随时间变化的特性曲线转化为幅值随频率变化的特性曲线
系统设计：利用傅里叶变换的快速傅里叶变换特性（fft）
系统测试：可以得出信号函数的功率谱函数，从时、频两域分析信号
遇到的问题及解决方法：mag=abs(y)由于傅里叶变换是复数域的变换，需要对其函数进行求模
结论：利用傅里叶变换及其逆变换可以简单的将信号函数进行时、频域的转化，有利于进行信号分析。
本人有一定从事信号分析的经历，并且积累了一定的经验，对傅里叶变换有比较深的了解，希望对你有所帮助，如果有什么问题可以在线问我。
欢迎分享，转载请注明来源：内存溢出
微信扫一扫
支付宝扫一扫
评论列表（0条）