xinhaochuli.rar_信号功率谱图_信号处理

clear all; clc load('yy15') lag=100; fs=100; y=yy15([1:160]); %选取前160点 %画出所给数据前160点的图 figure(1) plot(y,'-') title('数据前160点图'); xlabel('t/Ts'); ylabel('x(t)'); [c,lags]=xcorr(y,lag,'unbiased'); Pxx=fft(c,length(lags)); %利用FFT变换计算信号的功率谱 fp=(0:length(Pxx)-1)'*fs/length(Pxx); %求功率谱的横坐标f Pxmag=abs(Pxx); %求幅值 figure(2) plot(fp(1:length(Pxx)/2),Pxmag(1:length(Pxx)/2)); %绘制功率谱曲线 title('功率谱图'); xlabel('f/Hz');ylabel('|Pxx|'); %画出原信号及频谱图 N=length(yy15(1:4096)); t=(0:N-1); figure(3) subplot(2,1,1);plot(yy15);title('时域图形'); y=fft(yy15); f=(0:1/N:1/2-1/N)*fs; subplot(2,1,2);plot(f,abs(y(1:N/2)));grid;xlabel('hz');title('频谱'); %第一个频率10Hz f2=8;f3=12;f4=14;M=1000; wc1=2*f2/fs;wc2=2*f3/fs;wc3=2*f4/fs; %归一化角频率，用于设置通带与阻带范围 f1=[0 wc1-0.05 wc1 wc2 wc2+0.05 1]; A=[0 0 1 1 0 0]; %设置带通或带阻，1为带通，0为带阻 weigh=[1 1 1 ]; %设置通带和阻带的权重 b=remez(60,f1,A,weigh); %传函分子，滤波器阶数为60 h1=freqz(b,1,M); %幅频特性 figure(4) f=(0:1/M:1-1/M)*fs/2; subplot(2,1,1);plot(f,abs(h1));grid;title('10Hz带通滤波器'); x1=filter(b,1,yy15); S1=fft(x1); f=(0:1/N:1/2-1/N)*fs; subplot(2,1,2);plot(f,abs(S1(1:N/2)));grid;xlabel('f/hz');title('滤波器频谱特性'); figure(5) plot(x1(1:160)); title('10Hz信号前160点时域图') grid %第二个频率22Hz f21=21;f22=22;f23=23; wc21=2*f21/fs;wc22=2*f23/fs;wc23=2*f23/fs; f2=[0 wc21-0.06 wc21 wc22 wc22+0.06 1]; A=[0 0 1 1 0 0]; weigh=[1 1 1 ]; b2=remez(60,f2,A,weigh); h21=freqz(b2,1,M); figure(6) f=(0:1/M:1-1/M)*fs/2; subplot(2,1,1);plot(f,abs(h21));grid;title('22Hz带通滤波器'); x2=filter(b2,1,yy15); S2=fft(x2); f=(0:1/N:1/2-1/N)*fs; subplot(2,1,2);plot(f,abs(S2(1:N/2)));grid;xlabel('f/hz');title('滤波器频谱特性') figure(7) plot(x2(1:160)); title('22Hz信号前160点时域图') grid %第三个频率30Hz f31=29.5;f32=30;f33=30.5; wc31=2*f31/fs;wc32=2*f32/fs;wc33=2*f33/fs; f3=[0 wc31-0.06 wc31 wc32 wc32+0.06 1]; A=[0 0 1 1 0 0]; weigh=[1 1 1 ]; b3=remez(60,f3,A,weigh); h31=freqz(b3,1,M); figure(8) f=(0:1/M:1-1/M)*fs/2; subplot(2,1,1);plot(f,abs(h31));grid;title('30Hz带通滤波器'); x3=filter(b3,1,yy15); S3=fft(x3); f=(0:1/N:1/2-1/N)*fs; subplot(2,1,2);plot(f,abs(S3(1:N/2)));grid;xlabel('f/hz');title('滤波器频谱特性') figure(9) plot(x3(1:160)); title('30Hz信号前160点时域图') grid %还原原信号 x=x1+x2+x3; figure(10); subplot(2,1,1); plot(x(1:160)); title('叠加后信号前160点') xlabel('t/Ts'); ylabel('x(t)'); grid on subplot(2,1,2); plot(yy15(1:160)); title('原始信号前160点') xlabel('t/Ts'); ylabel('x(t)'); grid on