Demonstration of Non-parametric Power Spectrum Estimation

The graph shows spectrum of data passed through a filter calculated using pwelch MATLAB command to evaluate the periodogram by welches method.


%Design nice filter
Wp = [.2 .7], Ws = [.1 .8];Rp=1,Rs=30;

[N,Wn]=ellipord(Wp,Ws,Rp,Rs);

[b,a]=ellip(N,Rp,Rs,Wn);

N=2^12;n=0:N-1;

% generate and filter data

%x=rand(size(n))*2-1;

x=randn(size(n));

y=filter(b,a,x);

%Periodogram

[Px,f] = pwelch(x,rectwin(N),0,N,1);

[Py,f] = pwelch(y,rectwin(N),0,N,1);

plot(f,10*log10(Px),f,10*log10(Py))

hold on

title('Periodogram rectangular window');xlabel('frequency cycles per sample');ylabel('P')

% Bartlett Rectangular window

[Px,f] = pwelch(x,rectwin(N/16),0,N/16,1);

[Py,f] = pwelch(y,rectwin(N/16),0,N/16,1);

plot(f,10*log10(Py),'g','LineWidth',3)

% Bartlett Hamming window

[Px,f] = pwelch(x,hamming(N/16),0,N/16,1);

[Py,f] = pwelch(y,hamming(N/16),0,N/16,1);

plot(f,10*log10(Py),'c','LineWidth',2)

legend('x:periodogram','y:periodogram','y:Bartlett,rect','y:Bartlett,hamming')

xlabel('frequency cycles per sample');ylabel('P')

[H,w]=freqz(b,a);plot(w/2/pi,20*log10(abs(H)),'r')

axis([0 .5 -80 20]);grid