welch.m

给出了用MATLAB语言实现功率谱估计的程序

%两种加窗分段平均周期图法
function y=welch()
%求信号x(n)=sin(2*pi*f1*n)+2cos(2*pi*f2*n)+w(n), n=1,···,N-1的功率谱。其中f1=0.05,f2=0.12,w(n)为白噪声
%N=1024
clf;
N=1024;
Ns=256;
n=[0:N-1];
f1=0.05;
f2=0.12;
wn=randn(1,N)
xn=sin(2*pi*f1*n)+2*cos(2*pi*f2*n)+wn;
%无重叠分段、加窗平均周期图法
w=hanning(256)';
Pxx1=abs(fft(w.*xn(1:256),Ns).^2)/norm(w)^2;
Pxx2=abs(fft(w.*xn(257:512),Ns).^2)/norm(w)^2;
Pxx3=abs(fft(w.*xn(513:768),Ns).^2)/norm(w)^2;
Pxx4=abs(fft(w.*xn(769:1024),Ns).^2)/norm(w)^2;
Pxx=10*log10((Pxx1+Pxx2+Pxx3+Pxx4)/4);
f=[0:length(Pxx)-1]/length(Pxx);
figure(3);
plot(f,Pxx);
xlabel('Frequency');
ylabel('Power Spectrum(dB)');
title('Averaged Modified Periodogram (no overlap) N=4*256');
grid;
%==========================================
% 2:1 重叠分段、加窗平均周期图法
w=hanning(256)';

Pxx1=abs(fft(w.*xn(1:256),Ns).^2)/Ns;
Pxx2=abs(fft(w.*xn(129:384),Ns).^2)/norm(w)^2;
Pxx3=abs(fft(w.*xn(257:512),Ns).^2)/norm(w)^2;;
Pxx4=abs(fft(w.*xn(385:640),Ns).^2)/norm(w)^2;;
Pxx5=abs(fft(w.*xn(513:768),Ns).^2)/norm(w)^2;;
Pxx6=abs(fft(w.*xn(641:896),Ns).^2)/norm(w)^2;;
Pxx7=abs(fft(w.*xn(769:1024),Ns).^2)/norm(w)^2;;
Pxx=10*log10((Pxx1+Pxx2+Pxx3+Pxx4+Pxx5+Pxx6+Pxx7)/7);
f=[0:length(Pxx)-1]/length(Pxx);
figure(4);
plot(f,Pxx);
xlabel('Frequency');
ylabel('Power Spectrum(dB)');
title('Averaged Modified Periodogram (halfoverlap) N=1024');
grid;