psd_881.m

来自「用dsp解压mp3程序的算法」· M 代码 · 共 32 行

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% psd_811.m - Power spectrum estimation program used in Section 8.1.1
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

load in.dat -ascii;  % load the signal in in.dat
NFFT = 2^(nextpow2(length(in))); % find the next power of 2 length
NFFT1 = 256;	% window of 256-sample is used for averaged periodogram
Fs = 10000;
figure;

% PSD plot (1): Periodogram

[Pxx1,F1] = periodogram(in,[],[],Fs);
subplot(3,1,1); plot(F1,10*log10(Pxx1)); 
title('Power spectrum of sine waves (Periodogram)');
grid on;axis([0 5000 -100 50])

% PSD plot (2): Bartlett using rectangular window and no overlapping

[Pxx2,F2] = pwelch(in,boxcar(NFFT1),0,NFFT1,Fs);
subplot(3,1,2); plot(F2,10*log10(Pxx2)); 
title('Power spectrum of sine waves (Bartlett)');
grid on; axis([0 5000 -100 50])

% PSD plot (3): Welch using Hamming window and 50% overlapping

[Pxx3,F3] = pwelch(in,hamming(NFFT1),NFFT1/2,NFFT1,Fs);
subplot(3,1,3); plot(F3,10*log10(Pxx3));
title('Power spectrum of sine waves (Welch)');
grid on; axis([0 5000 -100 50])