samp9_7.m

MATLAB7.x数字信号处理 光盘内容

%Samp9_7
clf;Fs=1000;   %采样频率
N=1024;Nsec=256;n=0:N-1;t=n/Fs; %数据长度，分段数据长度、时间序列
w=hanning(256)';   %采用的窗口数据
randn('state',0);   %设置产生随机数的初始状态
xn=sin(2*pi*50*t)+2*sin(2*pi*120*t)+randn(1,N);  %带噪声的信号
%采用不重叠加窗方法的功率谱估计
pxx1=abs(fft(w.*xn(1:256),Nsec).^2)/norm(w)^2; %第一段加窗振幅谱平方
pxx2=abs(fft(w.*xn(257:512),Nsec).^2)/norm(w)^2; %第二段加窗振幅谱平方
pxx3=abs(fft(w.*xn(513:768),Nsec).^2)/norm(w)^2; %第三段加窗振幅谱平方
pxx4=abs(fft(w.*xn(769:1024),Nsec).^2)/norm(w)^2; %第四段加窗振幅谱平方
Pxx=10*log10((pxx1+pxx2+pxx3+pxx4)/4);    %求得平均功率谱，转换为dB
f=(0:length(Pxx)-1)*Fs/length(Pxx);  %求得频率序列
subplot(2,1,1),plot(f(1:Nsec/2),Pxx(1:Nsec/2)); %绘制功率谱曲线
xlabel('频率/Hz');ylabel('功率谱/dB');
title('加窗平均周期图(无重叠) N=4*256');
grid on
%采用重叠加窗方法的功率谱估计
pxx1=abs(fft(w.*xn(1:256),Nsec).^2)/norm(w)^2; %第一段加窗振幅谱平方
pxx2=abs(fft(w.*xn(129:384),Nsec).^2)/norm(w)^2; %第二段加窗振幅谱平方
pxx3=abs(fft(w.*xn(257:512),Nsec).^2)/norm(w)^2; %第三段加窗振幅谱平方
pxx4=abs(fft(w.*xn(385:640),Nsec).^2)/norm(w)^2; %第四段加窗振幅谱平方
pxx5=abs(fft(w.*xn(513:768),Nsec).^2)/norm(w)^2; %第五段加窗振幅谱平方
pxx6=abs(fft(w.*xn(641:896),Nsec).^2)/norm(w)^2; %第六段加窗振幅谱平方
pxx7=abs(fft(w.*xn(769:1024),Nsec).^2)/norm(w)^2; %第七段加窗振幅谱平方
Pxx=10*log10((pxx1+pxx2+pxx3+pxx4+pxx5+pxx6+pxx7)/7);%平均功率谱转换为dB
f=(0:length(Pxx)-1)*Fs/length(Pxx); %频率序列
subplot(2,1,2),plot(f(1:Nsec/2),Pxx(1:Nsec/2)); %绘制功率谱曲线
xlabel('频率/Hz');ylabel('功率谱/dB');
title('加窗平均周期图(重叠一半)N=1024');
grid on