数据采集与信号处理.docx
《数据采集与信号处理.docx》由会员分享，可在线阅读，更多相关《数据采集与信号处理.docx（21页珍藏版）》请在冰豆网上搜索。
数据采集与信号处理
哈尔滨理工大学
研究生考试试卷
考试科目：
数据采集与信号处理
阅卷人：
专业：
姓名：
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计算。
处理结果为采用窗口显示时域波形和频域波形。
此信号的时域谱，频域谱，功率谱如下图所示：
其MATLAB代码为：
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）信号经过均值化处理或不经过均值化处理的结果比较
通过以上两个图的分析，我们可以看出均值化处理后的频谱的低频段消失，这就去去除了常规的干扰频谱，如环境噪声等，对我们进行频谱分析有很大作用。
其MATLAB代码为
%含直流分量而未均值化的信号
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）采用不同窗函数时的谱结果（矩形窗函数,汉宁窗函数，汉明窗）
其MATLAB代码为：
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）典型函数的频谱（矩形窗函数,汉宁窗函数，直线，阶跃函数，δ函数，方波，三角波等）
此部分MATLAB代码如下：
t=0:
0.001:
0.2;
N=256;
FS=300;
w=boxcar（N）;%产生信号
subplot（211）;plot（w）;title（'矩形窗函数的时域波形图'）;
axis（[0,260,0,2]）;gridon;
y=fft（w,N）;%FFT运算
mag=abs（y）;%取幅值
f=（0:
length（y）-1）*FS/length（y）;
subplot（212）;plot（f（1:
N/2）,mag（1:
N/2））;%输出FS/2点幅频谱图
title（'矩形窗函数频域波形图'）;grid;
xlabel（'频率'）;ylabel（'幅值'）;
t=0:
0.001:
0.2;
N=256;
FS=300;
w=hanning（N）;%产生信号
subplot（211）;plot（w）;title（'汉宁窗函数的时域波形图'）;gridon;
y=fft（w,N）;%FFT运算
mag=abs（y）;%取幅值
f=（0:
length（y）-1）*FS/length（y）;
subplot（212）;
plot（f（1:
N/2）,mag（1:
N/2））;%输出FS/2点幅频谱图
title（'汉宁窗函数频域波形图'）;grid;
xlabel（'频率'）;ylabel（'幅值'）;
t=0:
0.001:
0.2;
N=256;
FS=300;
w=1;%产生信号
y=fft（w,N）;%FFT运算
mag=abs（y）;%取幅值
f=（0:
length（y）-1）*FS/length（y）;
plot（f（1:
N/2）,mag（1:
N/2））;%输出FS/2点幅频谱图
title（'直线频域波形图'）;grid;
xlabel（'Frequency（Hz）'）;ylabel（'Magnitude'）;
%阶跃函数的频域波形图
clc;clf;t=0:
0.001:
0.2;
N=256;
FS=300;
w=ones（1,N）;%产生信号
subplot（211）;plot（w）;
title（'阶跃函数的时域波形图'）;grid;
y=fft（w,N）;%FFT运算
mag=abs（y）;%取幅值
f=（0:
length（y）-1）*FS/length（y）;
subplot（212）;
plot（f（1:
N/2）,mag（1:
N/2））;%输出FS/2点幅频谱图
title（'阶跃函数的频域波形图'）;grid;
xlabel（'频率'）;ylabel（'幅值'）;
t=0:
0.001:
0.2;
N=256;
FS=300;
w=zeros（1,N）;w
（1）=1;%产生信号
subplot（211）;plot（w）;grid;
title（'δ函数的时域波形图'）;
y=fft（w,N）;
升级会员 