我设计的是多相插值滤波器,算法自己推的。
clear;
[n,Wn,beta,ftype]=kaiserord([1/8 5/24],[1 0],[0.00008 0.00008],1); %得出kaiser窗的参数n=62,beta=8.0708
fs=1;fc=1/6;w=2*fc/fs;%将模拟滤波器的技术指标转换为数字滤波器的技术指标.由滤波器的参数要求得出所需滤波器的阶数为n+1=63阶,符合3倍插值3的整数倍的需要,故无须再定义n值
window=kaiser(n+1,beta); %使用kaiser窗函数
b=fir1(n,[w],window); %用标准频率响应的加窗设计函数fir1,并得到滤波器系数
freqz(b); %数字滤波器频率响应
由此得滤波器的系数。
clear;
h=[0.0000 -0.0000 -0.0001 -0.0002 0.0000 0.0004 0.0005 -0.0000 -0.0010 -0.0014 0.0000 0.0024 0.0030 -0.0000 -0.0047 -0.0059 0.0000 0.0088 0.0106 -0.0000 -0.0153 -0.0184 0.0000 0.0267 0.0324 -0.0000 -0.0499 -0.0647 0.0000 0.1357 0.2746 0.3333 0.2746 0.1357 0.0000 -0.0647 -0.0499 -0.0000 0.0324 0.0267 0.0000 -0.0184 -0.0153 -0.0000 0.0106 0.0088 0.0000 -0.0059 -0.0047 -0.0000 0.0030 0.0024 0.0000 -0.0014 -0.0010 -0.0000 0.0005 0.0004 0.0000 -0.0002 -0.0001 -0.0000 0.0000];
%设置滤波器系数
[n,Wn,beta,ftype]=kaiserord([75/200 77/200],[1 0],[0.0001 0.0001],1);
fs=1;fc=76/200;w=2*fc/fs;
window=kaiser(n+1,beta); %使用kaiser窗函数
b=fir1(n,[w],window); %用标准频率响应的加窗设计函数fir1,并得到滤波器系数
t=1:1:1024;
s=sin(256*t)+sin(512*t)+sin(1024*t);
x=filter(b,1,s);
w=zeros(3072,1);
n=1;
m=1;
while m<3072
w(m)=x(n);
m=m+1;
w(m)=0;
m=m+1;
w(m)=0;
m=m+1;
n=n+1;
end%输入x经插值后得到w
y1=zeros(1024,1);
y2=zeros(1024,1);
y3=zeros(1024,1);
for m=1:1:1024
if m<=21
for k=1:1:m
y1(m)=y1(m)+x(m-k+1)*h(3*k-2);
y2(m)=y2(m)+x(m-k+1)*h(3*k-1);
y3(m)=y3(m)+x(m-k+1)*h(3*k);
end
else
for k=1:1:21
y1(m)=y1(m)+x(m-k+1)*h(3*k-2);
y2(m)=y2(m)+x(m-k+1)*h(3*k-1);
y3(m)=y3(m)+x(m-k+1)*h(3*k);
end
end
end%输入x经过多相运算后得到三个输出y1,y2,y3
n=1;
m=1;
while n<3072
y(n)=y1(m);
n=n+1;
y(n)=y2(m);
n=n+1;
y(n)=y3(m);
n=n+1;
m=m+1;
end%将y1,y2,y3整合成输出y
fs=3072;N=3072;
S=abs(fft(y));
S = S./max(S);
S = 20*log10(S+1e-15);
xvals = linspace(0,floor(fs/2)-floor((fs/N)),N/2);
plot(xvals(1:N/2),S(1:N/2)), grid on,axis([0,floor(fs/2)-1,-150,20]);
xlabel('Frequency (Hz)');
ylabel('Magnitude (dB)');
衰减达不到80dB,我试了好几种方法,包括改变输入,改变滤波器阶数,都不行。
那位大侠帮帮我!