%清除命令窗口及变量
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]);
%作printf怎么实现的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;
———————————————————————————————————————
这是一个傅里叶变化的子函数，你可以自己做主函数传递你这里的参数验证
// 入口参数：
// 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;
}
FFT in C code 用C语言实现FFT
#include <stdio.h>
#include <math.h>
//此代码来源《数字信号处理C语言程序集》殷福亮、宋爱军，沈阳：辽宁科学技术出版社，1997.7
//数组x存储时域序列的实部，数组y存储时域序列的虚部
//n代表N点FFT，sign=1为FFT，sign=-1为IFFT
void FFT(double x[],double y[],int n,int sign)
{
int i,j,k,l,m,n1,n2;
double c,c1,e,s,s1,t,tr,ti;
//Calculate i = log2N
for(j = 1,i = 1; i<16; i++)
{
m = i;
j = 2*j;
if(j == n)
break;
}
//计算蝶形图的输入下标（码位倒读）
n1 = n - 1;
for(j=0,i=0; i<n1; i++)
{
if(i<j)
{
tr = x[j];
ti = y[j];
x[j] = x[i];
y[j] = y[i];
x[i] = tr;
y[i] = ti;
}
k = n/2;
while(k<(j+1))
{
j = j - k;
k = k/2;
}
j = j + k;
}
//计算每一级的输出，l为某一级，i为同一级的不同，使用同一内存（即位运算）
版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。
发表评论