clc;clear;close all;
t=0:0.0001:0.5;
f1=2000;f2=1000;f3=4000;f4=5600;
xt=sin(f1*2*pi*t)+sin(f2*2*pi*t)+sin(f3*2*pi*t)+sin(f4*2*pi*t);  
plot(1000*t(1:100),xt(1:100));title('抽樣前模擬信號(hào)');
grid;
pause;

%抽樣
Fs=45000;  %抽樣頻率
L=301;
n=0:L-1;
xn=sin(f1/Fs*2*pi*n)+sin(f2/Fs*2*pi*n)+sin(f3/Fs*2*pi*n)+sin(f4/Fs*2*pi*n);
xnM=max(xn);
axis([0,L-1,-1.5*xnM,1.5*xnM]);
bar(n,xn,0);title('抽樣后的數(shù)字信號(hào)');grid;
pause;

%加窗后計(jì)算FFT
N=1024;   %FFT  點(diǎn)數(shù)
L=161;     %信號(hào)點(diǎn)數(shù)
x(1:L)=xn(31:(30+L));    %截?cái)嗉哟?x((L+1):N)=0;                  %補(bǔ)零
Xk=abs(fft(x));
XkM=max(Xk);
axis([0 10 0 1.1]);
f=(0:(N/2-1))*Fs/N;
% 換算為模擬頻率
Xkh=Xk(1:(N/2))/XkM;
plot(f,Xkh);grid;
title('輸入模擬信號(hào)x(t)的頻譜分析');
pause;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


fpass=4000;fstop=5000;fs=20000;
Apass=0.3; Astop=20;
wpass=2*pi*fpass;wstop=2*pi*fstop;
W=0:5:2*pi*fs/2;
leng=length(W);     %模擬原型濾波器指標(biāo)  
Wpass=fpass*2*pi/fs; Wstop=fstop*2*pi/fs;   %數(shù)字濾波器指標(biāo)
omigapass=tan(Wpass/2);omigastop=tan(Wstop/2);  
Epass=sqrt(10^(Apass/10)-1);Estop=sqrt(10^(Astop/10)-1);
e=Estop/Epass;w=omigastop/omigapass;     %雙線性變換后的模擬等效濾波器指標(biāo)
N=ceil(log(e+sqrt(e^2-1))/log(w+sqrt(w^2-1))); %N=6
K=floor(N/2);
a=1/N*log(1/Epass+sqrt(1+1/(Epass^2)));
for i=1:leng
    Hm(i)=sqrt(1/(1+Epass^2*Cn_cheb1(N,(W(i)/wpass))^2));   %公式11.6.46
end  

plot(W/(2*pi*1000),abs(Hm));
ylabel('|H(f)|');
xlabel('freq/kHz');
title('Type 1 chebyshev1,N=6');
grid;
pause;

%傳輸函數(shù)
for i=1:K     %thita
    q(i)=pi/(2*N)*(N-1+2*i);%公式11.16.17
end   

for i=1:K
    s(i)=omigapass*sinh(a)*cos(q(i))+j*omigapass*cosh(a)*sin(q(i));%公式11.6.53
end

if mod(N,2)==0
    s0=1;   %因?yàn)镹為even公式11.6.54
else
    s0=-omigapass*sinh(a);
end

omiga0=omigapass*sinh(a);
for i=1:K
    omiga(i)=omigapass*sin(q(i));
end      
if mod(N,2)==0
    Hs0=sqrt(1/(1+Epass^2));
else
   Hs0=omiga0/(s+omiga0);
end
omiga=0:0.01:pi;
lengomiga=length(omiga);
Hs=ones(lengomiga)*Hs0;
for k=1:lengomiga
    for i=1:K
    Hs(k)=Hs(k)*(omiga0^2+omiga(i)^2)/(-omiga(k)^2-2*omiga0*cos(q(i))*j*omiga(k)+omiga0^2+omiga(i)^2);
    end
end

plot(omiga/pi,abs(Hs));%雙線性變換的等效模擬濾波器
title('模擬等效濾波器幅頻響應(yīng)');
grid on;
pause;
%計(jì)算數(shù)字濾波器參數(shù)
%寫出系統(tǒng)函數(shù)
%G0=0.1790
%G=[ 0.0409    0.0351    0.0326]
%a1=[-1.6417   -1.4048   -1.2962]
%a2=[  0.8055    0.5453    0.4267]
%H0=0.9441
for i=1:K;
    G(i)=(omiga0^2+omiga(i)^2)/(1-2*omiga0*cos(q(i))+omiga0^2+omiga(i)^2);
    a1(i)=2*(omiga0^2+omiga(i)^2-1)/(1-2*omiga0*cos(q(i))+omiga0^2+omiga(i)^2);
    a2(i)=(1+2*omiga0*cos(q(i))+omiga0^2+omiga(i)^2)/(1-2*omiga0*cos(q(i))+omiga0^2+omiga(i)^2);
end
G0=omiga0/(omiga0+1);
a01=(omiga0-1)/(omiga0+1);%公式11.6.65-11.6.67

if mod(N,2)==0
    Hz0=sqrt(1/(1+Epass^2));
%else
%    Hz0=(G0*(1+1/z))/(1+a01/z);
end
w=0:0.05:pi;
z=exp(j*w);
ll=length(z);
Hz=ones(ll)*Hz0;
for k=1:ll
    for i=1:K
    Hz(k)=Hz(k)*G(i)*(1+1/z(k))^2/(1+a1(i)/z(k)+a2(i)/z(k)^2);
    end
end
plot(w/pi,abs(Hz));%數(shù)字濾波器的頻響 
title('數(shù)字濾波器頻率響應(yīng)');
grid on;
pause;


%離散信號(hào)通過數(shù)字濾波器
Lengsignal=1001;
b0=G;b1=2*G;b2=G;
a01=a1(1);a11=a1(2);a21=a1(3);
a02=a2(1);a12=a2(2);a22=a2(3);


b00=b0(1);b10=b0(2);b20=b0(3);
b01=b1(1);b11=b1(2);b21=b1(3);
b02=b2(1);b12=b2(2);b22=b2(3);

%狀態(tài)初始化
w00=0;w10=0;w20=0;
w01=0;w11=0;w21=0;
w02=0;w12=0;w22=0;

%輸入信號(hào)
for i=0:Lengsignal-1
    x=sin(f1/Fs*2*pi*i)+sin(f2/Fs*2*pi*i)+sin(f3/Fs*2*pi*i)+sin(f4/Fs*2*pi*i)
    x0=x;
    w00=-a01*w01-a02*w02+x0;
    y0=b00*w00+b01*w01+b02*w02;
    w02=w01;w01=w00;
    
    x1=y0;
    w10=x1-a11*w11-a12*w12;
    y1=b10*w10+b11*w11+b12*w12;
    w12=w11;w11=w10; 
    
    x2=y1;
    w20=x2-a21*w21-a22*w22;
    y(i+1)=b20*w20+b21*w21+b22*w22;
    w22=w21;w21=w20;
end 

%輸出序列
n=0:Lengsignal-1;
yM=max(y);
axis([0 L-1 -1.1*yM 1.1*yM]);
bar(100*n,y(1:Lengsignal),0)
xlabel('n')
ylabel('y(n)')
title('輸出濾波后的數(shù)字信號(hào)')
pause;


NO=1024;
lengout=201;
yn(1:lengout)=y(21:(20+lengout));
yn((lengout+1):NO)=0;  %補(bǔ)零
Yk=abs(fft(yn));
YkM=max(Yk);
axis([0 10 -40 1]);
f=(0:(NO/2-1))*Fs/NO;
Ykh=Yk(1:(NO/2))/YkM;
plot(f,20*log10(Ykh));
xlabel('f in kHz');
ylabel('Y(2*pi*f) in dB');
title('輸出序列頻譜分析/dB');
grid;
pause;
plot(f,Ykh)%輸出序列的頻域分析
title('輸出序列頻譜分析');
grid;
pause;
t1=0:1000;
plot(t1(1:400)/Fs,y(1:400));%輸出的模擬信號(hào)
title('輸出模擬信號(hào)');
grid;