功率谱估计实验.doc

Power Spectrum Estimation上机实验 姓名：戚永前 学号：10720885一：设输入音频信号，取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所示：取样fs=20KHzXa(t)2:1分段M=32三角窗加权规范运算M/2 点ODFT取模平方再除MU补奇数谱线内存UBxx(k)L个分段求平均D(2k)图2 2:1分段覆盖测量流程图计算机编程流程图如图3所示。

2:1分原信号为L段，M=32，三角窗加权每段，并将其重组为矩阵xn2逐段规范运算取实部、虚部，用ODFT得偶数谱线并求模平方除N由ODFT谱线共轭奇偶对称性补奇数谱线Pxx^(k)输入N及Ts开始结束图3 计算机测量流程图实验结果如图4所示,具体源程序参见附录程序源代码2图4 2:1分段覆盖所得的功率谱由图我们可以看到分段所得的功率谱分辨率较低，这与理论结果是一致的附录：程序源代码1(matlab)：%filename power_spectrum_estimation.m%version:10282009 pfx%this matlab program is used for power spectrum estimation. clear;clf;Fs=40000;%sample frequency 40kHz;N=128;%the number of samplesFc=2000;%the frequency of audio signal 2kHz;n=0:N-1;t=n/Fs;%time domainXt=cos(2*pi*Fc*t);%the audio signalTRGt=triang(N);%triangular windowU=sum(TRGt.^2)/N;% the energy of the windowYt=Xt.*TRGt';% the windowed signalfigure(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);endDk1=fft(Ut);%N/2 points fft ODFTDk2=fft(Vt);%N/2 points fft ODFTPs1=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 componentsendpxx1=10*log10(pow1);%display as dBpxx2=10*log10(pow2);%display as dBf=n/N*Fs;%the coordinate frequency subplot(2,2,3);plot(f,pxx1,'k');%display the power spectrum of the original signaltitle('Original Signal Power Spectrum Pxx(k) N=128');xlabel('frequency/Hz');ylabel('Power Spectrum(dB)');subplot(2,2,4);plot(f,pxx2,'k');%display the power spectrum of the windowed signaltitle('Windowed Signal Power Spectrum Bxx(k) N=128');xlabel('frequency/Hz');ylabel('Power Spectrum(dB)');程序源代码2(matlab)：%filename: power_spectrum_estimation_2nd.m%version: 10302009 pfx%this matlab program is used for power spectrum estimation. clear;clf;Fs=40000;%sample frequency 40kHz;N=128;%the number of samplesM=32;L=7;Fc=2000;%the frequency of audio signal 2kHz;n=0:N-1;t=n/Fs;% time domainXt=cos(2*pi*Fc*t);%the audio signalTRGt=triang(M);%triangular windowU=sum(TRGt.^2)/M;% the energy of the windowXt1=Xt(1:32);Xt2=Xt(17:48);Xt3=Xt(33:64);Xt4=Xt(49:80);Xt5=Xt(65:96);Xt6=Xt(81:112);Xt7=Xt(97:128);Yt1=Xt1.*TRGt';% the windowed signal1Yt2=Xt2.*TRGt';% the windowed signal2Yt3=Xt3.*TRGt';% the windowed signal3Yt4=Xt4.*TRGt';% the windowed signal4Yt5=Xt5.*TRGt';% the windowed signal5Yt6=Xt6.*TRGt';% the windowed signal6Yt7=Xt7.*TRGt';% the windowed signal7for m=0:M/2-1%normalization ODFT Ut1(m+1)=(Yt1(m+1)-j*Yt1(M/2+m+1))*exp(-j*m*pi/M); Ut2(m+1)=(Yt2(m+1)-j*Yt2(M/2+m+1))*exp(-j*m*pi/M); Ut3(m+1)=(Yt3(m+1)-j*Yt3(M/2+m+1))*exp(-j*m*pi/M); Ut4(m+1)=(Yt4(m+1)-j*Yt4(M/2+m+1))*exp(-j*m*pi/M); Ut5(m+1)=(Yt5(m+1)-j*Yt5(M/2+m+1))*exp(-j*m*pi/M); Ut6(m+1)=(Yt6(m+1)-j*Yt6(M/2+m+1))*exp(-j*m*pi/M); Ut7(m+1)=(Yt7(m+1)-j*Yt7(M/2+m+1))*exp(-j*m*pi/M);endDk1=fft(Ut1);%M/2 points fft ODFTDk2=fft(Ut2);%M/2 points fft ODFTDk3=fft(Ut3);%M/2 points fft ODFTDk4=fft(Ut4);%M/2 points fft ODFTDk5=fft(Ut5);%M/2 points fft ODFTDk6=fft(Ut6);%M/2 points fft ODFTDk7=fft(Ut7);%M/2 points fft ODFTPs1=abs(Dk1).^2/M/U;;%the power of the magnitude Ps2=abs(Dk2).^2/M/U;;%the power of the magnitude Ps3=abs(Dk3).^2/M/U;;%the power of the magnitude Ps4=abs(Dk4).^2/M/U;;%the power of the magnitude Ps5=abs(Dk5).^2/M/U;;%the power of the magnitude Ps6=abs(Dk6).^2/M/U;;%the power of the magnitude Ps7=abs(Dk7).^2/M/U;;%the power of the magnitude for m=0:M/2-1 pow1(2*m+1)=Ps1(m+1);%the even components pow1(M-2*m)=Ps1(m+1);%add the odd components pow2(2*m+1)=Ps2(m+1);%the even components pow2(M-2*m)=Ps2(m+1);%add the odd components pow3(2*m+1)=Ps3(m+1);%the even components pow3(M-2*m)=Ps3(m+1);%add the odd components pow4(2*m+1)=Ps4(m+1);%the even components pow4(M-2*m)=Ps4(m+1);%add the odd components pow5(2*m+1)=Ps5(m+1);%the even components pow5(M-2*m)=P。