哈尔滨理工大学
研究生考试试卷
考试科目：数据采集与信号处理阅卷人：
专业：
姓名：
2013年06月21日
一、基本内容：基于FFT的功率谱分析程序设计与应用
1．基本要求
1）对一个人为产生的信号进行采用FFT变换方法进行功率谱分析。
已知信号x(n)=80.0*COS(2*3.14*SF*n/FS)
式中：n=0,1,2 ……N-1
SF---信号频率
FS---采样频率
其FFT变换结果X(k)可用下面提供的FFT子程序求出，计算功率谱的公式为：
W（k）=2(XR(k)2 +XI(k)2)/N
式中：k=0,1,2 ……N/2-1
XR(k)--- X(k)的实部
XI(k)--- X(k)的虚部
请用VB,VC或C++Builder编译器编程，或采用MATLAB计算，或采用高级语言调用MATLAB计算。处理结果为采用窗口显示时域波形和频域波形。
此信号的时域谱，频域谱，功率谱如下图所示：
其MA TLAB代码为：
FS=200;
SF=10;
N=1024;
n=0:N-1;
t=n/FS;
x=80.0*cos(2*3.14*SF*t);
subplot(221);
plot(t,x);
xlabel('t');
ylabel('y');
title('x=80.0*cos(2*3.14*SF*t)时域波形');
grid;
y=fft(x,N);
mag=abs(y);
f=(0:length(y)-1)*FS/length(y);%进行对应的频率转换
subplot(222);
plot(f(1:N/2),mag(1:N/2));%做频谱图
xlabel('频率(Hz)');
ylabel('幅值');
title('x=80.0*cos(2*3.14*SF*t)幅频谱图N=1024');
grid;
Py =2*(y.*conj(y))/N; %计算功率谱密度Py
subplot(223)
plot(f(1:N/2),Py(1:N/2));
xlabel('频率(Hz)');
ylabel('功率谱密度');
title('x=80.0*cos(2*3.14*sf*t)功率谱密度');
grid;
二. 对一个用A/D数据采集板采集的信号进行频谱分析1）方波的频谱分析图像和程序
%fangbopufenxi
fid = fopen('F:\研究生信号处理\fanbo_45HZ_1024_1000HZ\fanbo _45HZ_1024_1000HZ '); %读入方波信号
SF=1000; %采样频率为1000HZ
[a,N]= fscanf(fid,'%f');
fclose(fid); %关闭打开的方波文件
y=fft(a,N); %进行快速傅里叶运算
Pyy =sqrt(y.*conj(y))*2.0/N; %取功率普密度
f=(0:length(Pyy)-1)*SF/length(Pyy);
LPyy=20*log10(Pyy);
plot(f(1:N/2),Pyy(1:N/2),'black'); %输出FS/2点幅频谱图grid;
2）三角波的频谱分析图像和程序
%sanjiaobopufenxi
fid = fopen('F:\研究生信号处理
\fanbo_45HZ_1024_1000HZ\sanjiao
_45HZ_1024_1000HZ'); %读入三角波信号
SF=1000; %采样频率为1000HZ
[a,N]= fscanf(fid,'%f');
fclose(fid); %关闭打开的三角波文件
y=fft(a,N); %进行快速傅里叶运算
Pyy =sqrt(y.*conj(y))*2.0/N; %取功率普密度
f=(0:length(Pyy)-1)*SF/length(Pyy);
LPyy=20*log10(Pyy);
plot(f(1:N/2),Pyy(1:N/2),'black'); %输出FS/2点幅频谱图grid;
3）正弦波的频谱分析图像和程序
%zhengxianbopufenxi
fid = fopen('F:\研究生信号处理
\fanbo_45HZ_1024_1000HZ\sin_45HZ_1024_1000HZ'); %读入三角波信号
SF=1000; %采样频率为1000HZ
[a,N]= fscanf(fid,'%f');
fclose(fid); %关闭打开的三角波文件
y=fft(a,N); %进行快速傅里叶运算
Pyy =sqrt(y.*conj(y))*2.0/N; %取功率普密度
f=(0:length(Pyy)-1)*SF/length(Pyy);
LPyy=20*log10(Pyy);
plot(f(1:N/2),Pyy(1:N/2),'black'); %输出FS/2点幅频谱图grid;
由上面的三幅图我们可以看出，在三角波、正弦波和方波这三种相同频率波的频谱分析中，方波的谐波特性最明显而且都是频率的奇数倍，三角波的谐波特性次之，正弦波的最不明显。
三、讨论
1) 信号经过均值化处理或不经过均值化处理的结果比较
20
4060801000
2000
4000
6000含有直流分量的余弦曲线未均值化的功率谱波形图频率(Hz)
幅值
20
406080
100
01000
20003000
4000含有直流分量的余弦曲线均值化后的功率谱波形图
频率(Hz)
幅值
通过以上两个图的分析，我们可以看出均值化处理后的频谱的低频段消失，这就去去除了常规的干扰频谱，如环境噪声等，对我们进行频谱分析有很大作用。
其MA TLAB 代码为
%含直流分量而未均值化的信号 FS=200; %采样频率 n=0:1:200;
SF=10; %信号频率
xn=80.0*cos(2*pi*SF*n/FS)+70; %产生波形序列 window=boxcar(length(xn)); %矩形窗 nfft=512; %采样点数
[Pxx,f]=periodogram(xn,window,nfft,FS); %直接法 subplot(211);
plot(f,Pxx);title('含有直流分量的余弦曲线未均值化的功率谱波形图'); xlabel('频率(Hz)'); ylabel('幅值'); grid;
%含直流分量而经均值化处理的信号 z=mean(xn); h=xn-z;
[Pxx1,f]=periodogram(h,window,nfft,FS); %直接法 subplot(212);
plot(f,Pxx1);title('含有直流分量的余弦曲线均值化后的功率谱波形图'); xlabel('频率(Hz)'); ylabel('幅值'); grid;
2）采用不同窗函数时的谱结果（矩形窗函数, 汉宁窗函数，汉明窗）
t
y
汉宁窗时域波形
20
406080
100
0100200300400500600
700800
900频率(Hz)
幅值
汉宁窗频域特性
0.2
0.4
0.6
0.8
1
1.2
1.4
-80-60-40-200204060
80t y
汉明窗时域波形
20
406080
100
050010001500200025003000
35004000
4500频率(Hz)
幅值
汉明窗频域特性
0.2
0.4
0.6
0.8
1
1.2
1.4
-80-60-40-200204060
80t y
矩形窗时域波形
20
406080
100
00.511.522.53
3.54
4.5x 10
4
频率(Hz)
幅值
矩形窗频域特性
其MA TLAB 代码为：
y=fft(x,N); mag=abs(y);
f=(0:length(y)-1)*FS/length(y);%进行对应的频率转换
w_han=(hanning(N))'; y1=x.*w_han; figure; plot(t,y1); xlabel('t'); ylabel('y');
title('汉宁窗时域波形'); grid;
y2=mag.*w_han; figure;
plot(f(1:N/2),y2(1:N/2)); xlabel('频率(Hz)'); ylabel('幅值');
title('汉宁窗频域特性'); grid;
w_rect=(rectwin(N))';
y3=x.*w_rect;
figure;
plot(t,y3);
xlabel('t');
ylabel('y');
title('矩形窗时域波形');
grid;
y4=mag.*w_rect;
figure;
plot(f(1:N/2),y4(1:N/2));
xlabel('频率(Hz)');
ylabel('幅值');
title('矩形窗频域特性');
grid;
w_ham=(hamming(N))';
y5=x.*w_ham;
figure;
plot(t,y5);
xlabel('t');
ylabel('y');
title('汉明窗时域波形');
grid;
y6=mag.*w_ham;
figure;
plot(f(1:N/2),y6(1:N/2));
xlabel('频率(Hz)');
ylabel('幅值');
title('汉明窗频域特性');
grid;
3）典型函数的频谱（矩形窗函数, 汉宁窗函数，直线，阶跃函数，δ函数，方波，三角波等）