完整word版C语言Matlab实现FFT几种编程实例Word文档格式.docx
- 文档编号:22195609
- 上传时间:2023-02-03
- 格式:DOCX
- 页数:14
- 大小:30.34KB
完整word版C语言Matlab实现FFT几种编程实例Word文档格式.docx
《完整word版C语言Matlab实现FFT几种编程实例Word文档格式.docx》由会员分享,可在线阅读,更多相关《完整word版C语言Matlab实现FFT几种编程实例Word文档格式.docx(14页珍藏版)》请在冰豆网上搜索。
nm1=FFT_N-1;
for(i=0;
i<
nm1;
i++)
{
if(i<
j)//如果i<
j,即进行变址
t=xin[j];
xin[j]=xin[i];
xin[i]=t;
}
k=nv2;
//求j的下一个倒位序
while(k<
=j)//如果k<
=j,表示j的最高位为1
{
j=j-k;
//把最高位变成0
k=k/2;
//k/2,比较次高位,依次类推,逐个比较,直到某个位为0
j=j+k;
//把0改为1
intle,lei,ip;
//FFT运算核,使用蝶形运算完成FFT运算
f=FFT_N;
for(l=1;
(f=f/2)!
=1;
l++)//计算l的值,即计算蝶形级数
;
for(m=1;
m<
=l;
m++)//控制蝶形结级数
{//m表示第m级蝶形,l为蝶形级总数l=log
(2)N
le=2<
<
(m-1);
//le蝶形结距离,即第m级蝶形的蝶形结相距le点
lei=le/2;
//同一蝶形结中参加运算的两点的距离
u.real=1.0;
//u为蝶形结运算系数,初始值为1
u.imag=0.0;
w.real=cos(PI/lei);
//w为系数商,即当前系数与前一个系数的商
w.imag=-sin(PI/lei);
for(j=0;
j<
=lei-1;
j++)//控制计算不同种蝶形结,即计算系数不同的蝶形结
for(i=j;
=FFT_N-1;
i=i+le)//控制同一蝶形结运算,即计算系数相同蝶形结
ip=i+lei;
//i,ip分别表示参加蝶形运算的两个节点
t=EE(xin[ip],u);
//蝶形运算,详见公式
xin[ip].real=xin[i].real-t.real;
xin[ip].imag=xin[i].imag-t.imag;
xin[i].real=xin[i].real+t.real;
xin[i].imag=xin[i].imag+t.imag;
u=EE(u,w);
//改变系数,进行下一个蝶形运算
/************************************************************
voidmain()
测试FFT变换,演示函数使用方法
无
************************************************************/
{
inti;
FFT_N;
i++)//给结构体赋值
s[i].real=sin(2*3.141592653589793*i/FFT_N);
//实部为正弦波FFT_N点采样,赋值为1
s[i].imag=0;
//虚部为0
FFT(s);
//进行快速福利叶变换
i++)//求变换后结果的模值,存入复数的实部部分
s[i].real=sqrt(s[i].real*s[i].real+s[i].imag*s[i].imag);
while
(1);
%////二、
%/////////////////////////////////////////////////////////////////////////////////////////////////////////////
%///////////////////////////////////////////////////////////////////////////////////////////////////////////
%////////////////////////////////MATLAB仿真信号源的源程序:
:
Clear;
Clc;
t=O:
O.01:
3;
yl=100*sin(pi/3*t);
n=l;
fort=-O:
10;
y2(1,n)=-61.1614*exp(-0.9*t);
n=n+;
end
min(y2)
y=[yl,y2];
figure
(1);
plot(y);
%funboxing(O.001+1/3)
%////////////////////////
%/////////快速傅里叶变换matlab程序:
%////////////////////////clc;
clear;
clf;
N=input('
Nodenumber'
)
T=input('
caiyangjiange'
f=input('
frenquency'
choise=input('
addzeroornot?
1/0'
n=0:
T:
(N-1)*T;
%采样点
k=0:
N-1;
x=sin(2*pi*f*n);
ifchoise==1
e=input('
Inputthenumberofzeros!
'
x=[xzeros(1,e)]
N=N+e;
else
end%加零
%给k重新赋值,因为有可能出现加零状况
bianzhi=bi2de(fliplr(de2bi(k,length(de2bi(N-1)))))+1;
%利用库函数进行变址运算
forl=1:
N
X(l)=x(bianzhi(l));
%将采样后的值按照变址运算后的顺序放入X矩阵中
d1=1;
form=1:
log2(N)
d2=d1;
%做蝶形运算的两个数之间的距离
d1=d1*2;
%同一级之下蝶形结之间的距离
W=1;
%蝶形运算系数的初始值
dw=exp(-j*pi/d2);
%蝶形运算系数的增加量
fort=1:
d2%
fori=t:
d1:
i1=i+d2;
if(i1>
N)break;
%判断是否超出范围
else
p=X(i1)*W;
X(i1)=X(i)-p;
X(i)=X(i)+p;
%蝶形运算
end
W=W*dw;
%蝶形运算系数的变化
endend
subplot(2,2,1);
t=0:
0.0000001:
N*T;
plot(t,sin(2*pi*f*t));
%画原曲线
subplot(2,2,2);
stem(k,x);
%画采样后的离散信号
subplot(2,2,3);
stem(k,abs(X)/max(abs(X)));
%画自己的fft的运算结果
subplot(2,2,4);
stem(k,abs(fft(x))/max(abs(fft(x))));
%调用系统的函数绘制结果,与自己的结果进行比较
//////三、
///////////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////FFT的C语言算法实现
////////////程序如下:
/************FFT***********/
#include<
stdio.h>
stdlib.h>
#defineN1000
typedefstruct
doublereal;
doubleimg;
}complex;
voidfft();
/*快速傅里叶变换*/
voidifft();
/*快速傅里叶逆变换*/
voidinitW();
voidchange();
voidadd(complex,complex,complex*);
/*复数加法*/
voidmul(complex,complex,complex*);
/*复数乘法*/
voidsub(complex,complex,complex*);
/*复数减法*/
voiddivi(complex,complex,complex*);
/*复数除法*/
voidoutput();
/*输出结果*/
complexx[N],*W;
/*输出序列的值*/
intsize_x=0;
/*输入序列的长度,只限2的N次方*/
doublePI;
intmain()
inti,method;
system("
cls"
);
PI=atan
(1)*4;
printf("
Pleaseinputthesizeofx:
\n"
/*输入序列的长度*/
scanf("
%d"
&
size_x);
Pleaseinputthedatainx[N]:
(suchas:
56)\n"
/*输入序列对应的值*/
size_x;
i++)
%lf%lf"
x[i].real,&
x[i].img);
initW();
/*选择FFT或逆FFT运算*/
UseFFT(0)orIFFT
(1)?
method);
if(method==0)
fft();
else
ifft();
output();
return0;
}
/*进行基-2FFT运算*/
voidfft()
inti=0,j=0,k=0,l=0;
complexup,down,product;
change();
log(size_x)/log
(2);
l=1<
i;
j+=2*l)
for(k=0;
k<
l;
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;
voidifft()
inti=0,j=0,k=0,l=size_x;
complexup,down;
(int)(log(size_x)/log
(2));
i++)/*蝶形运算*/
l/=2;
j+=2*l)
k++)
add(x[j+k],x[j+k+l],&
up.real/=2;
up.img/=2;
sub(x[j+k],x[j+k+l],&
down.real/=2;
down.img/=2;
divi(down,W[size_x*k/2/l],&
voidinitW()
W=(complex*)malloc(sizeof(complex)*size_x);
W[i].real=cos(2*PI/size_x*i);
W[i].img=-1*sin(2*PI/size_x*i);
voidchange()
complextemp;
unsignedshorti=0,j=0,k=0;
doublet;
k=i;
j=0;
t=(log(size_x)/log
(2));
while((t--)>
0)
j=j<
1;
j|=(k&
1);
k=k>
>
if(j>
i)
temp=x[i];
x[i]=x[j];
x[j]=temp;
voidoutput()/*输出结果*/
Theresultareasfollows\n"
%.4f"
x[i].real);
if(x[i].img>
=0.0001)
+%.4fj\n"
x[i].img);
elseif(fabs(x[i].img)<
0.0001)
%.4fj\n"
voidadd(complexa,complexb,complex*c)
c->
real=a.real+b.real;
img=a.img+b.img;
voidmul(complexa,complexb,complex*c)
real=a.real*b.real-a.img*b.img;
img=a.real*b.img+a.img*b.real;
voidsub(complexa,complexb,complex*c)
real=a.real-b.real;
img=a.img-b.img;
voiddivi(complexa,complexb,complex*c)
real=(a.real*b.real+a.img*b.img)/(
b.real*b.real+b.img*b.img);
img=(a.img*b.real-a.real*b.img)/(b.real*b.real+b.img*b.img);
/////四、
///////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
%//////选带傅里叶变换(zoom-fft)MATLAB
%移频(将选带的中心频率移动到零频)
%数字低通滤波器
(防止频率混叠)
%重新采样
(将采样的数据再次间隔采样,间隔的数据取决于分析的带宽,就是放大倍数)
%复FFT(由于经过了移频,所以数据不是实数了)
%频率调整(将负半轴的频率成分移到正半轴)
function
[f,
y]
=
zfft(x,
fi,
fa,
fs)
%
x为采集的数据
fi为分析的起始频率
fa为分析的截止频率
fs为采集数据的采样频率
f为输出的频率序列
y为输出的幅值序列(实数)
f0
(fi
+
fa)
/
2;
%中心频率
N
length(x);
%数据长度
r
0:
b
2*pi*f0.*r
./
fs;
x1
x
.*
exp(-1j
b);
%移频
bw
fa
-
fi;
B
fir1(32,
fs);
%滤波
截止频率为0.5bw
x2
filter(B,
1,
x1);
c
x2(1:
floor(fs/bw):
N);
N1
length(c);
f
linspace(fi,
N1);
y
abs(fft(c))
*
circshift(y,
[0,
floor(N1/2)]);
%将负半轴的幅值移过来
end
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////////////////
%上边四程序测试应用实例:
fs
2048;
T
100;
t
1/fs:
T;
30
cos(2*pi*110.*t)
cos(2*pi*111.45.*t)
25*cos(2*pi*112.3*t)
48*cos(2*pi*113.8.*t)+50*cos(2*pi*114.5.*t);
109,
115,
plot(f,
y);
图片效果:
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 完整 word 语言 Matlab 实现 FFT 编程 实例
![提示](https://static.bdocx.com/images/bang_tan.gif)