新建文件夹.zip_Piecewise Smooth_分段平滑_分段平滑的MATLAB代码实现_周期图平滑_谱线平滑

a=load('C:\Users\chenhao\Desktop\jianmo\2-1-4.txt'); Fs=4000; %采样频率 N=4096; %采样点数 y=mean(a); %两段平滑; Nsec=2048; %数据的长度和FFT所用的数据长度 Pxx1=abs(fft(a(1:2048)-y,Nsec).^2)/Nsec; %第一段功率谱 Pxx2=abs(fft(a(2049:4096)-y,Nsec).^2)/Nsec; %第二段功率谱 Pxx=10*log10((Pxx1+Pxx2)/2); %Fourier振幅谱平方的平均值，并转化为dB x1=std((Pxx1+Pxx2)/2)^2 f=(0:length(Pxx)-1)*Fs/length(Pxx); %给出频率序列 figure; plot(f(1:Nsec/2),Pxx(1:Nsec/2)) %绘制功率谱曲线 xlabel('频率/Hz'); ylabel('功率谱/dB'); title('N=2*2048'); grid on; %四段平滑; Nsec=1024; %数据的长度和FFT所用的数据长度 Pxx1=abs(fft(a(1:1024)-y,Nsec).^2)/Nsec; %第一段功率谱 Pxx2=abs(fft(a(1025:2048)-y,Nsec).^2)/Nsec; %第二段功率谱 Pxx3=abs(fft(a(2049:3072)-y,Nsec).^2)/Nsec; %第三段功率谱 Pxx4=abs(fft(a(3073:4096)-y,Nsec).^2)/Nsec; %第四段功率谱 Pxx=10*log10((Pxx1+Pxx2+Pxx3+Pxx4)/4); %Fourier振幅谱平方的平均值，并转化为dB x2=std((Pxx1+Pxx2+Pxx3+Pxx4)/4)^2 f=(0:length(Pxx)-1)*Fs/length(Pxx); %给出频率序列 figure; plot(f(1:Nsec/2),Pxx(1:Nsec/2)) %绘制功率谱曲线 xlabel('频率/Hz'); ylabel('功率谱/dB'); title('N=4*1024'); grid on; %八段平滑; Nsec=512; %数据的长度和FFT所用的数据长度 Pxx1=abs(fft(a(1:512)-y,Nsec).^2)/Nsec; %第一段功率谱 Pxx2=abs(fft(a(513:1024)-y,Nsec).^2)/Nsec; %第二段功率谱 Pxx3=abs(fft(a(1025:1536)-y,Nsec).^2)/Nsec; %第三段功率谱 Pxx4=abs(fft(a(1537:2048)-y,Nsec).^2)/Nsec; %第四段功率谱 Pxx5=abs(fft(a(2049:2560)-y,Nsec).^2)/Nsec; %第五段功率谱 Pxx6=abs(fft(a(2561:3072)-y,Nsec).^2)/Nsec; %第六段功率谱 Pxx7=abs(fft(a(3073:3584)-y,Nsec).^2)/Nsec; %第七段功率谱 Pxx8=abs(fft(a(3585:4096)-y,Nsec).^2)/Nsec; %第八段功率谱 Pxx=10*log10((Pxx1+Pxx2+Pxx3+Pxx4+Pxx5+Pxx6+Pxx7+Pxx8)/8); %Fourier振幅谱平方的平均值，并转化为dB x3=std((Pxx1+Pxx2+Pxx3+Pxx4+Pxx5+Pxx6+Pxx7+Pxx8)/8)^2 f=(0:length(Pxx)-1)*Fs/length(Pxx); %给出频率序列 figure; plot(f(1:Nsec/2),Pxx(1:Nsec/2)); %绘制功率谱曲线 xlabel('频率/Hz'); ylabel('功率谱/dB'); title('N=8*512'); grid on;