低通滤波器

function num=lxcfir(fp,fs,Fs,s) %fp为通带截止频率，fs阻带截止频率,Fs=1/T %s=1 for lowpass filter,s=-1 for highpass filter %sample hxcfir(25,20,100,1) wp=2*pi*fp/Fs; ws=2*pi*fs/Fs; if s==1 kuan=ws-wp; end if s==-1 kuan=wp-ws; end M=ceil(6.6*pi/kuan)+1;%确定窗宽为M %ceil(X) 求出离X最近的整数值 n=[0:1:M-1];%表示0-（M-1）这M个数 wc=(ws+wp)/2;%临界频率 %wc=central w if s==1 hd=lixiang(wc,M,1);%调用lixiang(wc,M,1) 子函数 end if s==-1 if rem(M,2)==0%M除以2的余数，即M为偶数 M=M+1;%变为奇数 end hd=lixiang(wc,M,-1); end window=(hamming(M))';%选出海明窗函数，"'"表示转置矩阵 h=hd.*window;%.* 表示数组乘法 fenzid=h;%分子 fenmud=1;%分母 t=0:0.01:1;%采样时间0.01 %输入信号 s1=sin(2*pi*t*6); s2=sin(2*pi*t*33); s=s1+s2; temp1=0; temp2=0; y=zeros(1,length(s));%length(x)=max(size(x))，返回矢量s的长度， %zeros(1,3)=(0,0,0) %滤波器的频率响应 %FIR滤波器的输出列可以按卷积写成：y(n)= h(0)x(n) +h(1)x(n-1)+...+h(N-1)x(n-(N-1)) %因为FIR是线性相位的，其冲击响应具有中心对称特点，上式可写成：y(n)=h(0)[x(0)+x(n-(N-1))]+h(1)[x(n-1)+x( %-(N-2))]+...+h((N-1)/2)[x(N-1)/2+x(N+1)/2] for i=1:length(s) if i<length(fenzid)%fenzid=h=numerator for j=1:i temp1=temp1+fenzid(j)*s(i-j+1); end else for j=1:length(fenzid) temp1=temp1+fenzid(j)*s(i-j+1); end end y(i)=temp1-temp2; temp2=0; if i<length(fenmud)-1 for k=1:i temp2=temp2+fenmud(k+1)*y(i-k+1); end else for k=1:length(fenmud)-1 temp2=temp2+fenmud(k+1)*y(i-k+1); end end temp1=0; end figure(1) subplot(2,1,1) plot(s); subplot(2,1,2) plot(y); %子函数程序 function h=lixiang(wc,N,s) M=(N-1)/2;%确定冲击相应位移系数，N为窗宽 for k=-M:M if k==0 h(k+M+1)=(1-s)/2+s*wc/pi; else h(k+M+1)=s*sin(wc*k)/(pi*k); end end