基于MATLAB的快速傅立叶分析程序设计.docx
《基于MATLAB的快速傅立叶分析程序设计.docx》由会员分享，可在线阅读，更多相关《基于MATLAB的快速傅立叶分析程序设计.docx（14页珍藏版）》请在冰豆网上搜索。
基于MATLAB的快速傅立叶分析程序设计
1.已知信号数据
对一个人为产生的信号进行采用FFT变换方法进行功率谱分析。
已知信号
x（n）=120.0*COS（2*3.14*SF*n/FS）
式中：
n=0,1,2……N-1
SF---信号频率
FS---采样频率
这里，定义参数如下：
fs=200；%设定采样频率
N=512；
sf=10；%设定余弦信号频率
采样点=1024；
2．信号的时域波形和频域波形
3．信号经过均值化处理或不经过均值化处理的结果比较
4．采用不同窗函数时的谱结果（矩形窗函数,汉宁窗函数，汉明窗等等）
5．典型函数的频谱（矩形窗函数,汉宁窗函数，三角窗函数，切比雪夫窗）
6．整周期和非整周期比较
基于MATLAB的快速傅立叶分析程序设计
1.已知信号数据
对一个人为产生的信号进行采用FFT变换方法进行功率谱分析。
已知信号x（n）=120.0*COS（2*3.14*SF*n/FS）
式中：
n=0,1,2……N-1SF---信号频率FS---采样频率
这里，定义参数如下：
fs=200；%设定采样频率
N=512；
sf=10；%设定余弦信号频率
采样点=1024；
2．信号的时域波形和频域波形
图1余弦函数的时域波形图
图2余弦函数的幅频谱波形
图3余弦函数的功率谱密度波形
上面三图的程序如下所示：
fs=200;%设定采样频率
N=512;
n=0:
N-1;
t=n/fs;
sf=10;%设定正弦信号频率
%生成信号
x=120.0*cos（2*3.14*sf*t）;
figure;
plot（t,x）;%作余弦信号的时域波形
xlabel（'t'）;
ylabel（'y'）;
title（'x=120.0*cos（2*3.14*sf*t）时域波形'）;
grid;
%进行FFT变换并做频谱图
y=fft（x,N）;%进行fft变换
mag=abs（y）;%求幅值
f=（0:
length（y）-1）'*fs/length（y）;%进行对应的频率转换
figure;
plot（f,mag）;%做频谱图
xlabel（'频率（Hz）'）;
ylabel（'幅值'）;
title（'x=120.0*cos（2*3.14*sf*t）幅频谱图N=512'）;
grid;
Py=2*（y.*conj（y））/N;%计算功率谱密度Py
figure;
plot（f,Py）;
xlabel（'频率（Hz）'）;
ylabel（'功率谱密度'）;
title（'x=120.0*cos（2*3.14*sf*t）功率谱密度'）;
grid;
3．信号经过均值化处理或不经过均值化处理的结果比较
图4余弦函数均值化处理后的时频域波形
结论：
经过均值化处理的波形发生上下移动，但是频谱不变。
4．采用不同窗函数时的谱结果（矩形窗函数,汉宁窗函数，汉明窗等等）
1）采用矩形窗函数时信号的频谱分析结果
图
图5矩形窗函数处理信号时域波形
图6矩形窗函数处理信号频域波形
2）采用汉明窗函数时信号的频谱分析结果
图7汉明窗函数处理信号时域波形
图8汉明窗函数处理信号频域波形
3）采用汉宁窗函数时信号的频谱分析结果
图9汉宁窗函数处理信号时域波形
图10汉宁窗函数处理信号频域波形
4）采用三角窗函数时信号的频谱分析结果
图11三角窗函数处理信号时域波形
图12三角窗函数处理信号频域波形
5）采用Blackman窗函数时信号的频谱分析结果
图13Blackman窗函数处理信号时域波形
图14Blackman窗函数处理信号频域波形
在前面程序的基础上，继续输入加窗处理程序，可以得到上面所示的结果，其加窗处理程序如下：
w_han=（hanning（N））';
y1=x.*w_han;
figure;
plot（t,y1）;
y2=mag.*w_han;
figure;
plot（f,y2）;
w_box=（boxcar（N））';
y3=x.*w_box;
figure;
plot（t,y3）;
y4=mag.*w_box;
figure;
plot（f,y4）
w_ham=（hamming（N））';
y5=x.*w_ham;
figure;
plot（t,y5）;
y6=mag.*w_ham;
figure;
plot（f,y6）
w_tri=（triang（N））';
y7=x.*w_tri;
figure;
plot（t,y7）;
y8=mag.*w_tri;
figure;
plot（f,y8）
w_black=（blackman（N））';
y9=x.*w_black;
figure;
plot（t,y9）;
y10=mag.*w_black;
figure;
plot（f,y10）
5．典型函数的频谱（矩形窗函数,汉宁窗函数，三角窗函数，切比雪夫窗）
设计方法：
主要应用了MALTAB中的交互式图形用户界面以及直接编程来做信号的处理过程，其设计过程如下：
1）图形用户界面的启动：
在MATLABCOMMAND窗口下，键入sptool，会弹出一个SPTool窗口。
2）在进行处理之前，我们需要建立一个所要处理的信号公式，即已知信号x（n）=120.0*cos（2*3.14*SF*n/FS），MATLAB提供了编程的代码，其代码如下：
Fs=200;%设定采样频率
N=512;
n=0:
N-1;t=n/Fs;
sf=10;%设定余弦信号频率
x=120.0*cos（2*3.14*sf*t）;%生成信号
plot（t,x）
grid
savehdata.matxFs%把已知信号保存到了工作空间中，以备调用
这样程序运行结果会生成信号数据文件hdata.mat，存放信号x和采样频率的数据。
3）从SPTool窗口中的File菜单中选择Import命令，弹出ImporttoSPTool窗。
窗口如下所示：
图15信号处理ImporttoSPTool窗口
设置好后，点击OK。
4）这样就可以进行信号的时域分析了，只需要点击下图中Signals中所导入的已知信号就可以了。
图16MATLAB信号处理SPTool窗口
5）如上图所示的Spectra为信号的频谱分析工具栏，选好所导入的已知信号，点击Crate按钮，进入频谱分析的窗口。
窗口图片如下所示：
图17信号的频谱分析工具栏Spectra窗口
通过上面图形所示，我们可以从左边的工具栏上选择各种窗函数，并且可以随意定义Nfft和Nwind的大小，方法比较简单直观。
同样，我们还可以用编程的方法实现上面要求，其各种窗函数及程序代码如下所示：
1）矩形窗函数
图18矩形窗函数频谱图
2）汉宁窗函数
图19汉宁窗函数频谱图
3）三角窗函数
图20三角窗函数频谱图
4）切比雪夫窗函数
图21切比雪夫窗函数频谱图
典型窗函数的频谱程序如下所示：
生成一个长度为50的矩形窗，并观察其频率特性
n=50;
window=boxcar（n）;
[h,w]=freqz（window,1）;
subplot（2,1,1）
stem（window）;
subplot（2,1,2）;
plot（w/pi,20*log（abs（h）/abs（h
（1））））;
生成一个长度为60的汉宁窗，并观察其频率特性
n=60;
window=hanning（n）;
[h,w]=freqz（window,1）;
subplot（1,2,1）
stem（window）;
subplot（1,2,2）;
plot（w/pi,20*log（abs（h）/abs（h
（1））））;
生成一个长度为40的三角窗，并观察其频率特性
n=40;
window=triang（n）;
[h,w]=freqz（window,1）;
subplot（1,2,1）
stem（window）;
subplot（1,2,2）;
plot（w/pi,20*log（abs（h）/abs（h
（1））））;
生成一个长度为40的切比雪夫窗，并观察其频率特性
n=50;
r=50;
window=chebwin（n,r）;
[h,w]=freqz（window,1）;
subplot（1,2,1）
stem（window）;
subplot（1,2,2）
plot（w/pi,20*log（abs（h）/abs（h
（1））））;
6．整周期和非整周期比较
以高度为1的矩形窗w（t）在一个完整周期内截取采样序列，在时域以w（t）截取采样序列，在频域内等效于W（f）和X（f）的卷积。
在f＝±m/T0,处正好是X（f）*W（f）=W（f-f0）的零点处，故在频谱上仅在f0处有单一频线，无‘泄漏”现象。
因此，只要矩形窗宽包含序列的若干完整周期，就都能在f'=±m/T0处使采样序列的加窗谱不歪曲其真实谱。
当用矩形窗函数w（t）对周期信号进行非整周期截取，如上图所示产生明显的“泄漏现象”，在W（f）*X（f）的谱图中，在f0≠nΔf，谱线的位置nΔf不尽为频谱的零点重合，W（f）的旁瓣会在频率为±m/T0处产生，使序列截断后的频谱和真正的频谱不同，有“泄漏”。
因此，在对信号进行周期采样时，采样长度最好是信号周期长度的整数倍，以防采集的信号失真。
非整周期信号处理分析：
图22非整周期信号时域分析
图23非整周期信号频域分析
升级会员 