tmp606.m

信号分析与处理的一些源代码

%MATLAB PROGRAM 6-6
clf
Fs=1000;
%Estimate PSD by the nonoverlapping signal section
%=================================================
N=1024;
Nsec=256;
n=0:N-1;
t=n/Fs;
xn=sin(2*pi*50*t)+2*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 (no overlap) N=4*256')
grid

%Estimate PSD by the half overlapping signal section
%===================================================
Fs=1000;
N=1024;
Nsec=256;
n=0:N-1;
t=n/Fs;
xn=sin(2*pi*50*t)+2*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