功率谱估计实验.doc

Power Spectrum Estimation上机实验 姓名：戚永前 学号一：设输入音频信号，取f=2KHz，fs=40KHz，N=128，W(n)为三角窗序列。用计算机求出功率谱估值及不分段的。 实验中使用matlab语言进行编程，在MATLAB R2008b上调试通过。所得到的结果如图1，具体源代码参见附录程序源代码1。 图1 实验1结果图 上图左下图所示的是所对应到实际频率的功率谱，右下图是不分段的所对应到实际频率的功率谱。图中我们可以清楚的看到，在2kHz频率处由一个尖峰，这与我们测试的音频信号的频率正好是吻合的，证明实验结果是正确的。 二：设输入音频信号，取f=2KHz，fs=40KHz，N=128，W(n)为三角窗序列。用计算机求分段的。用2:1覆盖分段，设各段的长度为M=32。 具体测量流程图如图2所示： 图2 2:1分段覆盖测量流程图 计算机编程流程图如图3所示。 图3 计算机测量流程图 实验结果如图4所示,具体源程序参见附录程序源代码2。 图4 2:1分段覆盖所得的功率谱 由图我们可以看到分段所得的功率谱分辨率较低，这与理论结果是一致的。 附录： 程序源代码1(matlab)： %filename power_spectrum_estimation.m %versionpfx %this matlab program is used for power spectrum estimation. clear; clf; Fs=40000;%sample frequency 40kHz; N=128;%the number of samples Fc=2000;%the frequency of audio signal 2kHz; n=0:N-1; t=n/Fs;%time domain Xt=cos(2*pi*Fc*t);%the audio signal TRGt=triang(N);%triangular window U=sum(TRGt.^2)/N;% the energy of the window Yt=Xt.*TRGt;% the windowed signal figure(1); subplot(2,2,1); plot(t,Xt,k);%display the original audio signal; title(Original Audio signal , Fs=40kHz,N=128,Fc=2kHz); xlabel(t/s); ylabel(X(t)/V); subplot(2,2,2); plot(t,Yt,k);%display the windowed audio signal; title(Windowed Audio signal, N=128,Triangular Window); xlabel(t/s); ylabel(Y(t)/V); for m=0:N/2-1%normalization ODFT Ut(m+1)=(Xt(m+1)-j*Xt(N/2+m+1))*exp(-j*m*pi/N); Vt(m+1)=(Yt(m+1)-j*Yt(N/2+m+1))*exp(-j*m*pi/N); end Dk1=fft(Ut);%N/2 points fft ODFT Dk2=fft(Vt);%N/2 points fft ODFT Ps1=abs(Dk1).^2;%the power of the magnitude Ps1=Ps1/N; Ps2=abs(Dk2).^2;%the power of the magnitude Ps2=Ps2/N/U; for m=0:N/2-1 pow1(2*m+1)=Ps1(m+1);%the even components pow1(N-2*m)=Ps1(m+1);%add the odd components pow2(2*m+1)=Ps2(m+1);%the even components pow2(N-2*m)=Ps2(m+1);%add the odd components end pxx1=10*log10(pow1);%display as dB pxx2=10*log10(pow2);%display as dB f=n/N*Fs;%the coordinate frequency subplot(2,2,3); plot(f,pxx1,k);%display the power spectrum of the or