%Samp10_2
Fs=100;dt=1/Fs;   %給定模擬輸入信號采樣間隔
f1=6;f2=10;f3=19;   %模擬輸入信號的頻率成分
N=500;             %數據點數
t=[0:N-1]*dt;   %模擬輸入信號的時間序列
x=sin(2*pi*f1*t)+0.5*cos(2*pi*f2*t)+0.5*sin(2*pi*f3*t); %模擬輸入信號
X=fft(x);   %得到輸入信號的Fouier變換
%設計寬帶Butterworth模擬濾波器
wp=[1 40]*2*pi;ws=[0.2 42.5]*2*pi; %通帶和阻帶截止頻率
Rp=1;Rs=20;   %通帶波紋和阻帶衰減
[Order,Wn]=buttord(wp,ws,Rp,Rs,'s');   %求解Butterworth濾波器的最小階數
w=linspace(0,Fs/2,N/2)*2*pi;   %設置繪制幅頻特性的頻率
[b1,a1]=butter(Order,Wn,'s');     %設計Butterworth帶通濾波器
H=freqs(b1,a1,w);   %計算帶通濾波器的頻率響應
magH=abs(H);phaH=unwrap(angle(H));  %求得振幅譜和相位譜
figure(1)
subplot(3,1,1),plot(w/(2*pi),20*log10(magH));  %繪制幅頻特性
xlabel('頻率/Hz')
ylabel('振幅/dB');
ylim([-100 0]);xlim([0,50]); grid on
title('寬帶模擬帶通濾波器');

%模擬輸出
H=freqs(b1,a1,w);   %濾波器的幅頻響應
H1=zeros(1,N);
for ii=1:N/2   %構建與Fourier變換得到前后數據共軛的形式
    H1(ii)=H(ii);    
    H1(N-ii)=conj(H1(ii));
end
X=fft(x);   %將原信號進行Fourier變換
Y=zeros(1,N);
for ii=1:N
    Y(ii)=X(ii).*H1(ii);   %按（10-4）式得到的輸出
end
y=real(ifft(Y));     %將頻率域的數據通過Fourier變換轉換到時間域
subplot(3,1,2),plot(t,x),title('輸入信號')  %繪制輸入信號
subplot(3,1,3),plot(t,y)  %繪制輸出信號
title('寬帶濾波器的模擬輸出信號') 
xlabel('時間/s')

%恢復地面運動
XX=zeros(1,N);
for ii=1:N
    if(abs(H1(ii))>1.0e-1)   %為防止位于阻帶內數據不穩定
        XX(ii)=Y(ii)./H1(ii);   %恢復地面運動公式（10-7）
     end
end
xx=real(ifft(XX));  %轉換到時間域
figure(2)
plot(t,x,t,xx,':r')
legend('原始信號','恢復信號',1)  %給出圖例
title('原始信號和恢復信號')
xlabel('時間/s')
%窄帶濾波器的設計
wp=[5 7]*2*pi;ws=[3.5 7.5]*2*pi;  %窄帶濾波器的通帶和阻帶邊界頻率
Rp=1;Rs=20;   %窄帶濾波器的通帶波紋和阻帶衰減
[Order,Wn]=buttord(wp,ws,Rp,Rs,'s'); %計算窄帶濾波器的最小階數和截止頻率  
w=linspace(0,Fs/2,N/2)*2*pi;   %設置繪制窄帶濾波器幅頻特性的頻率
[b2,a2]=butter(Order,Wn,'s');   %設計Butterworth窄帶濾波器
H=freqs(b2,a2,w);    %計算Butterworth窄帶濾波器的頻率響應
magH=abs(H);phaH=unwrap(angle(H));  %計算振幅譜和相位譜
figure(3)
subplot(3,1,1),plot(w/(2*pi),20*log10(magH));  %繪制幅頻特性
xlabel('頻率/Hz'),ylabel('振幅/dB'); ylim([-100 0]);xlim([0,50])
title('窄帶模擬濾波器');
grid on
%模擬輸出
H=freqs(b2,a2,w);
H2=zeros(1,N);
for ii=1:N/2
    H2(ii)=H(ii);
    H2(N-ii)=conj(H2(ii));
end
X=fft(x);
Y1=zeros(1,N);
for ii=1:N
    Y1(ii)=X(ii).*H2(ii);
end
y1=real(ifft(Y1));
subplot(3,1,2),plot(t,x),title('輸入信號')  %繪制輸入信號
subplot(3,1,3),plot(t,y1)  %繪制輸出信號
title('窄帶濾波器的模擬輸出') 
xlabel('時間/s')
%運用寬頻帶儀器的輸出仿真得到窄帶儀器上的輸出
figure(4)
XY=zeros(1,N);
for ii=1:N
    if(abs(H1(ii))>1.0e-2)
        XY(ii)=Y(ii)*H2(ii)/H1(ii);  %公式（10-6）
    end
end
xy=real(ifft(XY));   %轉換到時間域
plot(t,y1,t,xy,':r')
legend('窄帶輸出','仿真輸出',1)
xlabel('時間/s')