clf;Fs=1000;%数据采样点数
%运用信号不重叠分段估计功率谱
N=1024;Nsec=256;n=0:N-1;t=n/Fs;%数据点数,分段时间,时间序列
randn('state',0);%设置产生随机数的初始状态
xn=sin(2*pi*100*t)+2*sin(2*pi*200*t)+3*sin(2*pi*400*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(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(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);%功率谱平均并转化为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
评论0