哈尔滨理工大学研究生考试试卷考试科目：数据采集与信号处理阅卷人：专业：姓名：2013年06月21日一、基本内容：基于FFT的功率谱分析程序设计与应用1．基本要求1）对一个人为产生的信号进行采用FFT变换方法进行功率谱分析。

已知信号x(n)=80.0*COS(2*3.14*SF*n/FS)式中：n=0,1,2 ……N-1SF---信号频率FS---采样频率其FFT变换结果X(k)可用下面提供的FFT子程序求出，计算功率谱的公式为：W（k）=2(XR(k)2 +XI(k)2)/N式中：k=0,1,2 ……N/2-1XR(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; %计算功率谱密度Pysubplot(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）方波的频谱分析图像和程序%fangbopufenxifid = 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）三角波的频谱分析图像和程序%sanjiaobopufenxifid = 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）正弦波的频谱分析图像和程序%zhengxianbopufenxifid = 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) 信号经过均值化处理或不经过均值化处理的结果比较204060801000200040006000含有直流分量的余弦曲线未均值化的功率谱波形图频率(Hz)幅值2040608010001000200030004000含有直流分量的余弦曲线均值化后的功率谱波形图频率(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）采用不同窗函数时的谱结果（矩形窗函数, 汉宁窗函数，汉明窗）00.20.40.60.81 1.2 1.4ty汉宁窗时域波形204060801000100200300400500600700800900频率(Hz)幅值汉宁窗频域特性0.20.40.60.811.21.4-80-60-40-20020406080t y汉明窗时域波形20406080100050010001500200025003000350040004500频率(Hz)幅值汉明窗频域特性0.20.40.60.811.21.4-80-60-40-20020406080t y矩形窗时域波形2040608010000.511.522.533.544.5x 104频率(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）典型函数的频谱（矩形窗函数, 汉宁窗函数，直线，阶跃函数，δ函数，方波，三角波等）此部分MA TLAB代码如下：t=0:0.001:0.2;N=256;FS=300;w=boxcar(N); %产生信号subplot(211);plot(w);title('矩形窗函数的时域波形图');axis([0,260,0,2]);grid on;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('汉宁窗函数的时域波形图');grid on;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); %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=square(2*pi*50*t); %产生信号subplot(211);plot(t,w);title('方波的时域波形图');axis([0,0.2,-0.2,1.2]);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=sawtooth(2*pi*50*t,0.5); subplot(211);plot(t,w); grid;title('三角波的时域波形图');%产生信号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('幅值');4）整周期和非整周期采样时两者的比较00.20.40.60.81-100-50050100ty整周期时域波形0100200300400500600246x 104频率(Hz)幅值整周期)幅频谱图N=102400.20.40.60.81 1.2 1.4-100-50050100ty非整周期时域波形050100150200250300350400450244频率(Hz)幅值非整周期幅频谱图N=1024此部分的MATLAB 代码为： FS=900; SF=10; N=1024; n=0:N-1; t=n/FS;x=80*cos(2*3.14*SF*t);subplot(211); plot(t,x); xlabel('t'); ylabel('y');title('非整周期时域波形'); grid;y=fft(x,N); mag=abs(y);f=(0:length(y)-1)*FS/length(y);%进行对应的频率转换subplot(212);plot(f(1:N/2),mag(1:N/2));%做频谱图xlabel('频率(Hz)');ylabel('幅值');title('非整周期幅频谱图N=1024');grid;FS=1024;SF=10;N=1024;n=0:N-1;t=n/FS;x=80.0*cos(2*3.14*SF*t);subplot(211);plot(t,x);xlabel('t');ylabel('y');title('整周期时域波形');grid;y=fft(x,N);mag=abs(y);f=(0:length(y)-1)*FS/length(y);%进行对应的频率转换subplot(212);plot(f(1:N/2),mag(1:N/2));%做频谱图xlabel('频率(Hz)');ylabel('幅值');title('整周期)幅频谱图N=1024');grid;5）讨论实验结果通过以上分析，采用整周期采样能够去除频谱中的干扰成分，容易得到准确的频率。