%
% 算法工匠 信号源的产生和滤波1
% firdesign.m
% author:照着抄的 作者是up主。
% 2018年6月3日 18:12:35
clear ;
close all;
clc;

fc1 = 10;
fc2 = 100;
fc3 = 450;  %三个频率分量

fs = 1000;      %采样频率
point_s =  1000; %采样点数
time = [1:point_s]/fs;   %时域抽样时间变量
delta_f = 1*fs/point_s;    %频率分辨率
sin_s1 = 0.5*sin( 2*pi*fc1*time);
sin_s2 = 1.1*sin( 2*pi*fc2*time);
sin_s3 = 0.8*sin( 2*pi*fc3*time);
sin_s = sin_s1 + sin_s2 + sin_s3; %生成及合成信号
fft_s1 = fft(sin_s);
figure(1);
subplot(2,1,1);
plot(time,sin_s);
title(\'原始时域信号\');
f = 0:point_s/2-1;
f = f*delta_f;
danbianfudu = 2*abs(fft_s1(1:point_s/2 ));
subplot(2,1,2);
plot(f,danbianfudu);
title(\'频域幅度特性\');

%%%%%%滤波处理
%设计低通滤波器，其上限截止频率为fc = 120Hz
%  wc = fc/(fs/2);
N = 128;    %滤波器阶数
fc = 120;
wc = fc/(fs/2); 
fir1_lowpass_filter =  fir1(N,wc);
figure(2);
freqz(fir1_lowpass_filter); %查看幅度特性和相位特性(归一化频率)
title(\'128阶低通滤波器\');
%实现滤波
filter_sin_s1 = filter(fir1_lowpass_filter,1,sin_s);
fft_low_filter = fft(filter_sin_s1); 
figure(3);
subplot(2,1,1);
plot(time,filter_sin_s1);
title(\'低通滤波后时域波形\');
subplot(2,1,2);
danbianfudu = 2*abs(fft_low_filter(1:point_s/2 ));
plot(f,danbianfudu);
title(\'低通滤波后的频率幅度特性\');
%滤波处理
%设计一个高筒滤波器，截止频率为15Hz
fc = 20 ;
wc = fc/(fs/2); %视频课件中之所以写0.15，是因为这个例子wc=0.15 对应于75Hz也能滤掉75Hz的，我这里取20Hz
fir1_hp_filter = fir1(N,wc,\'high\');
figure(4);
freqz(fir1_hp_filter);
filter_sin_s2 = filter(fir1_hp_filter,1,sin_s);
fft_high_filter = fft(filter_sin_s2);
figure(5);
subplot(2,1,1);
plot(time,filter_sin_s2);
title(\'高通滤波后的时域波形\');
subplot(2,1,2);
danbianfudu = 2*abs(fft_high_filter(1:point_s/2 ));
plot(f,danbianfudu);
title(\'高通滤波后的频域幅度特性\'); 

%带通滤波器
%只让100Hz的通过
%wc1 = 50/(fs/2);
%wc2 = 400/(fs/2);
fc1 = 50;
fc2 = 400;
wc1 = fc1 /( fs / 2 );
wc2 = fc2 /( fs / 2 );
fir1_band_pass_filter = fir1(N,[wc1 wc2],\'bandpass\');  %ftype
%
figure(6);
freqz(fir1_band_pass_filter);
%实现滤波
fir1_bp_s1 = filter(fir1_band_pass_filter,1,sin_s);
fft_bp_filter = fft(fir1_bp_s1);
figure(7);
subplot(2,1,1);
plot(time ,fir1_bp_s1 );
title(\'带通滤波器滤波后波形\');
subplot(2,1,2);
danbianfudu = 2*abs(fft_bp_filter(1:point_s/2 ));
plot(f,danbianfudu);
title(\'带通滤波后的频率幅度特性\');


%带阻滤波器
%只不让100Hz的通过
%wc1 = 50/(fs/2);
%wc2 = 400/(fs/2);
fc1 = 50;
fc2 = 400;
wc1 = fc1 /( fs / 2 );
wc2 = fc2 /( fs / 2 );
fir1_band_stop_filter = fir1(N,[wc1 wc2],\'stop\');  %ftype
%
figure(8);
freqz(fir1_band_stop_filter);
%实现滤波
fir1_bs_s1 = filter(fir1_band_stop_filter,1,sin_s);
fft_bs_filter = fft(fir1_bs_s1);
figure(9);
subplot(2,1,1);
plot(time ,fir1_bs_s1 );
title(\'带阻滤波器滤波后波形\');
subplot(2,1,2);
danbianfudu = 2*abs(fft_bs_filter(1:point_s/2 ));
plot(f,danbianfudu);
title(\'带阻滤波后的频率幅度特性\');