yelv123 发表于 2007-6-22 14:05

数字滤波器matlab的源代码

function lvbo(Ua,Ub,choise)
%参考指令：lvbo(2*pi,10*pi,1/0/-1)
U1=min(Ua,Ub);
U2=max(Ua,Ub);
Us=16*U2;
T=2*pi/Us;
T_sum=4*max(2*pi/Ua,2*pi/Ub);
sum=T_sum/T;
t=T:T:T_sum;
x=sin(U1*t)+0.8*sin(U2*t);
X=DFT(x);
figure(1); subplot(221)
U=Us/sum:Us/sum:Us;
stem(U,abs(X));grid on
axis()
title('原模拟信号采样频谱图')
Ucd=U1+(U2-U1)*1/5;Usd=U2-(U2-U1)*1/5;
switch choise
case 1
Hz_ejw=IIR_DF_BW(Ucd,1,Usd,30,T,sum);
case -1
Hz_ejw=IIR_DF_CF(Ucd,1,Usd,30,T,sum);
case 0
Hz_ejw=FIR_DF_HM(U1,U2,T,sum);
otherwise
Hz_ejw=IIR_DF_BW(Ucd,1,Usd,30,T,sum);
end
Y=X.*Hz_ejw;
y=1/sum*conj(DFT(conj(Y)));
figure(1); subplot(224)
plot(t,real(y)); title('模拟信号滤波后');grid on
axis()
subplot(222);
plot(t,x); hold on
axis()
x=sin(U1*t);
plot(t,x,':r');grid on
title('模拟信号滤波前')
function Hz_ejw=IIR_DF_BW(Ucd,Ap,Usd,As,t,sum)
% 巴特沃思滤波器
E=(10^(0.1*Ap)-1)^0.5;
V=(10^(0.1*As)-1)^0.5;
Wc=Ucd*t; Ws=Usd*t;
Ucd=Wc/t; Usd=Ws/t;
Uca=(2/t)*tan(Ucd*t/2); Usa=(2/t)*tan(Usd*t/2);
N=ceil(log10(V/E)/log10(Usa/Uca));
k=;
Spk=exp(j*(pi/2+(2*k-1)/(2*N)*pi));
i=find(real(Spk)<0);
Sk(1:N)=Spk(i);
den=real(poly(Sk'));
k0=polyval(den,0);
disp('模拟巴特沃思滤波器的归一化统函数 Ha(s) 为')
tf(k0,den)
syms s z T;
den_jU=1;
s=s/Uca;
for i=1:N
den_jU=s^(N-i+1)*den(i)+den_jU;
end
Ha_s=simple(1/den_jU);
H_z=subs(Ha_s,'s',(2/T)*((1-1/z)/(1+1/z)));
k=1:sum;w=(2*pi/sum)*k;
ejw=exp(j*w);
Hz_ejw=subs(H_z,{z,T},{ejw,t*ones(1,length(ejw))});
figure(1); subplot(223)
plot(w,abs(Hz_ejw)); grid on
title('巴特沃思低通滤波器')
axis()
function Hz_ejw=IIR_DF_CF(Ucd,Ap,Usd,As,t,sum)
% 切比雪夫低通滤波器
E=(10^(0.1*Ap)-1)^0.5;
V=(10^(0.1*As)-1)^0.5;
Wc=Ucd*t; Ws=Usd*t;
Ucd=Wc/t; Usd=Ws/t;
Uca=(2/t)*tan(Ucd*t/2); Usa=(2/t)*tan(Usd*t/2);
N=ceil(acosh(V/E)/acosh(Usa/Uca));;
A=1/E+(1/E^2+1)^0.5;
a=1/2*(A^(1/N)-A^(-1/N));
b=1/2*(A^(1/N)+A^(-1/N));
k=1:2*N;
Spk=-a*sin((2*k-1)/(2*N)*pi)+j*b*...
cos((2*k-1)/(2*N)*pi);
i=find(real(Spk)<0);
Sk(1:N)=Spk(i);
den=real(poly(Sk'));
k0=1;
disp('模拟切比雪夫低通滤波器的归一化统函数 Ha(s) 为')
tf(k0,den)
if (rem(N,2)==1)
for i=1:N
k0=k0*(-Sk(i));
end
elseif ((rem(N,2))==0)
k0=1;
for i=1:N
k0=k0*(-Sk(i));
end
end
if (rem(N,2)==0)
k0=10^(-0.05*Ap)*k0;
end
k0=real(k0);
syms s z T;
den_jU=1;
s=s/Uca;
for i=1:N
den_jU=s^(N-i+1)*den(i)+den_jU;
end
Ha_s=simple(1/den_jU);
H_z=subs(Ha_s,'s',(2/T)*((1-1/z)/(1+1/z)));
k=1:sum;w=(2*pi/sum)*k;
ejw=exp(j*w);
Hz_ejw=subs(H_z,{z,T},{ejw,t*ones(1,length(ejw))});
figure(1); subplot(223)
plot(w,abs(Hz_ejw));grid on
title('切比雪夫低通滤波器')
axis()
function Hz_ejw=FIR_DF_HM(U1,U2,T,sum)
wp=U1*T;ws=U2*T;
kuan=ws-wp;
M=sum;
n=;
wc=(ws+wp)/2;
hd=H_D(wc,M);
window=hamming_m(M);
h_z=hd.*window;
Hz_ejw=DFT(h_z);
k=1:sum;w=(2*pi/sum)*k;
figure(1); subplot(223)
plot(w,abs(Hz_ejw));grid on
axis();
title('海明窗函数低通滤波器')
function hd=H_D(wc,N)
M=(N-1)/2;
for k=-M:M
if k==0
hd(k+M+1)=wc/pi;
else
hd(k+M+1)=sin(wc*k)/(pi*k);
end
end
function wn=hamming_m(M)
n=0:M-1;
wn(n+1)=0.54-0.46*cos((2*pi*n)/(M-1));
function Xk=DFT(xn)
% 离散傅立叶变换,xn为原序列,Xk为DFT变换后的序列
N=length(xn);
n=0:N-1;
k=0:N-1;
WN=exp(-j*2*pi/N);
nk=n'*k;
WNnk=WN.^nk;
Xk=xn*WNnk;

qiuqia17 发表于 2007-11-21 15:26

Undefined command/function 'tf'

求解函数
:@L

graduate 发表于 2007-11-26 01:02

很好，支持:victory:

ilovechina1985 发表于 2008-5-6 23:35

数字滤波器是不是fir 的比iir的要好？也就是说FIR用的比较多？麻烦哪位指点一下

haoyuenan 发表于 2009-2-23 16:04

iir滤波器需要的阶数比较多，但是效果比较好，受运行条件限制，一般没有那么长的时间做高阶的运算，所以都是用fir滤波器吧，不知道是不是这样，呵呵
页: [1]
查看完整版本: 数字滤波器matlab的源代码