声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 3791|回复: 4

[共享资源] 数字滤波器matlab的源代码

[复制链接]
发表于 2007-6-22 14:05 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
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([Us/sum,Us/2,0,1.2*max(abs(X))])
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([0,T_sum,-max(real(y))*1.5,max(real(y))*1.5])
subplot(222);
plot(t,x); hold on
axis([0,T_sum,-max(x)*1.2,max(x)*1.2])
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=[1:2*N];
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([2*pi/sum,pi,-0.2,1.2*max(abs(Hz_ejw))])
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([2*pi/sum,pi,-0.5,max(abs(Hz_ejw))])
function Hz_ejw=FIR_DF_HM(U1,U2,T,sum)
wp=U1*T;ws=U2*T;
kuan=ws-wp;
M=sum;
n=[0:1:M-1];
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([2*pi/sum,pi,-0.2,1.2*max(abs(Hz_ejw))]);
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;

评分

1

查看全部评分

回复
分享到:

使用道具 举报

发表于 2007-11-21 15:26 | 显示全部楼层
Undefined command/function 'tf'

求解  函数
:@L
发表于 2007-11-26 01:02 | 显示全部楼层
很好,支持:victory:
发表于 2008-5-6 23:35 | 显示全部楼层
数字滤波器是不是fir 的比iir的要好?也就是说FIR用的比较多?麻烦哪位指点一下
发表于 2009-2-23 16:04 | 显示全部楼层
iir滤波器需要的阶数比较多,但是效果比较好,受运行条件限制,一般没有那么长的时间做高阶的运算,所以都是用fir滤波器吧,不知道是不是这样,呵呵

评分

1

查看全部评分

您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-11 23:10 , Processed in 0.073078 second(s), 24 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表