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
评论0