平均周期图法.m

来自「matlab,功率谱估计」· M 代码 · 共 41 行

clf;
%4分段平均周期图法
Fs=1000;
N=1024;
Nsec=256;
n=0:N-1;
t=n/Fs;
xn=sin(2*pi*50*t)+sin(2*pi*120*t)+randn(1,N);
Pxx1=abs(fft(xn(1:256),Nsec).^2)/Nsec;
Pxx2=abs(fft(xn(257:512),Nsec).^2)/Nsec;
Pxx3=abs(fft(xn(513:768),Nsec).^2)/Nsec;
Pxx4=abs(fft(xn(769:1024),Nsec).^2)/Nsec;
Pxx=10*log10((Pxx1+Pxx2+Pxx3+Pxx4)/4);
f=(0:length(Pxx)-1)*Fs/length(Pxx);
subplot(211)
plot(f,Pxx);
xlabel('Frequency (Hz)');ylabel('Power Spectrum (dB)');
title('Averaged Periodogram (none overlap) N=4*256');
grid
%4分段（2：1重叠）平均周期图法
Fs=1000;
N=1024;
Nsec=256;
n=0:N-1;
t=n/Fs;
xn=sin(2*pi*50*t)+sin(2*pi*120*t)+randn(1,N);
Pxx1=abs(fft(xn(1:256),Nsec).^2)/Nsec;
Pxx2=abs(fft(xn(129:384),Nsec).^2)/Nsec;
Pxx3=abs(fft(xn(257:512),Nsec).^2)/Nsec;
Pxx4=abs(fft(xn(385:640),Nsec).^2)/Nsec;
Pxx5=abs(fft(xn(513:768),Nsec).^2)/Nsec;
Pxx6=abs(fft(xn(641:896),Nsec).^2)/Nsec;
Pxx7=abs(fft(xn(769:1024),Nsec).^2)/Nsec;
Pxx=10*log10((Pxx1+Pxx2+Pxx3+Pxx4+Pxx5+Pxx6+Pxx7)/7);
f=(0:length(Pxx)-1)*Fs/length(Pxx);
subplot(212)
plot(f,Pxx);
xlabel('Frequency (Hz)');ylabel('Power Spectrum (dB)');
title('Averaged Periodogram (half overlap) N=1024');
grid