plot(f,abs(X(1:(N+1)/2)));title('原信号频域图') xlabel('频率(Hz)'); ylabel('幅度'); figure(1);
subplot(4,1,2);
plot(t,y);title('加50hz工频信号时域图') xlabel('时间(s)'); ylabel('幅度'); figure(2);
subplot(4,1,2);
plot(f,abs(Y(1:(N+1)/2)));title('加50hz工频信号频域图') xlabel('频率(Hz)'); ylabel('幅度'); figure(1)
subplot(4,1,3);
plot(t,y1);title('加谐波后的时域图') xlabel('时间(s)'); ylabel('幅度'); figure(2)
subplot(4,1,3);
plot(f,abs(Y1(1:(N+1)/2)));title('加谐波后的频域图') xlabel('频率(Hz)'); ylabel('幅度'); figure(1);
subplot(4,1,4);
plot(t,y1);title('加白噪后的时域图'); figure(2);
subplot(4,1,4);
plot(f,abs(Y1(1:(N+1)/2)));title('加白噪后的频域图');
(3)滤除部分:
滤除50HZ信号后的时域图(blackman低通滤波器):
滤除50HZ信号后的频域图:
blackman低通滤波器滤波:
具体代码:
%原心电信号
xl=load('E:\\心电信号-新.txt'); x=xl(:,2); N=size(x,1); fs=1000;
t=(0:length(x)-1)/fs;
f=fs/N*(0:(N+1)/2-1)+1;%f=fs/N*(0:N-1); X=fft(x,N);
%加50hz工频干扰 y=x+sin(2*50*pi*t)'; Y=fft(y,N);
%blackman低通滤波器 wp=30/fs*pi;ws=40/fs*pi; B=ws-wp;
M=ceil(12*pi/B)-1;
bl=fir1(M,(ws+wp)/2/pi,'low',blackman(M+1)); yl=fftfilt(bl,y); Yl=fft(yl,N); figure(1);
subplot(3,1,1);
plot(t,x);title('原信号时域图') figure(2); subplot(3,1,1)
plot(f,abs(X(1:(N+1)/2)));%plot(f,abs(X)); title('原信号频域图')
figure(1);
subplot(3,1,2);
plot(t,y);title('加50hz工频信号时域图') figure(2);
subplot(3,1,2);
plot(f,abs(Y(1:(N+1)/2)));%plot(f,abs(Y)); title('加50hz工频信号频域图') figure(1);
subplot(3,1,3);
plot(t,yl);title('滤波后信号时域图') figure(2)
subplot(3,1,3);
plot(f,abs(Yl(1:(N+1)/2)));%plot(f,abs(Yl)); title('滤波后信号频域图') figure(3)
[h2,w2]=freqz(bl,1);
plot(w2/(2*pi)*fs,20*log10(abs(h2)));title('blackman低通滤波器频率响应图');
滤除50HZ信号后的时域图(切比雪夫滤波器):
滤除50HZ信号后的频域图(切比雪夫滤波器):