MATLAB的快速傅立叶分析程序设计

基于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.信号的时域波形和频域波形绘制
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;
在前面程序的基础上，继续输入加窗处理程序，可以得到上面所示的结果，其加窗处理程序如下：
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)
3.典型函数的频谱（矩形窗函数,汉宁窗函数，三角窗函数，切比雪夫窗）
设计方法：主要应用了MALTAB中的交互式图形用户界面以及直接编程来做信号的处理过程，其设计过程如下：
1）图形用户界面的启动：在MATLAB COMMAND窗口下，键入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
save hdata.mat x Fs %把已知信号保存到了工作空间中，以备调用
这样程序运行结果会生成信号数据文件hdata.mat，存放信号x和采样频率的数据。

3）从SPTool窗口中的File菜单中选择Import命令，弹出Import to SPTool窗。

设置好后，点击OK。

4）这样就可以进行信号的时域分析了，只需要点击下图中Signals中所导入的已知信号就可以了。

5）如上图所示的Spectra为信号的频谱分析工具栏，选好所导入的已知信号，点击Crate 按钮，进入频谱分析的窗口。

通过上面图形所示，我们可以从左边的工具栏上选择各种窗函数，并且可以随意定义Nfft和Nwind的大小，方法比较简单直观。

同样，我们还可以用编程的方法实现上面要求，其各种窗函数及程序代码如下所示：
典型窗函数的频谱程序如下所示：
生成一个长度为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))));。