matlab中的数字滤波没有给出预期的结果

Question

我正在尝试在 matlab 中设计和测试一个简单的数字高通滤波器。 我有两个脚本：第一个是设计过滤器，第二个是递归算法的实现

第一个脚本：过滤器的设计：

s=sym('s');
z=sym('z');
w=sym('w');
p=sym('p');
f=sym('f');
x=sym('x');
Pi=sym('Pi');
Ts=sym('Ts');
vc=sym('vc');

n=2;
Fspb=1/(s^2+sqrt(2)*s+1);   % simple low pass butterworth
Fs=subs(Fspb, s, 2*Pi*vc/s); % tranform to high pass butterworth with vc cutting freq
Fz=subs(Fs, s, 2/Ts*(z-1)/(z+1)); %bilinear transformation
Fp=subs(Fz, z, exp(Ts*p))
Fw=subs(Fp, p, 1i*w);
Ff=subs(Fw, w, 2*Pi*f);

vmax=5000;
vc_val=1000; %1000 Hz
vccont=1/Pi/Ts*tan(Pi*Ts*vc_val);
Ts_val=0.0001
fval=0:0.1:vmax;

pretty(expand(Fz))

[num, den]=numden(Fz);
cn=coeffs(num,z)
cd=coeffs(den,z)

它给了我系数（对于 z）和滤波器的频率响应：

cn =     
[ 1, -2, 1]

cd =
[ Pi^2*Ts^2*vc^2 - 2^(1/2)*Pi*Ts*vc + 1, 2*Pi^2*Ts^2*vc^2 - 2, Pi^2*Ts^2*vc^2 + 2^(1/2)*Pi*Ts*vc + 1] 

这是使用freqz函数显示的mfilter：

还有我的第二个脚本，在一个简单的正弦函数上实现过滤器和测试：

Pi=3.14;
v=500;
Ts=0.0001;
x=0:Ts:100000*Ts;

y=sin(x*2*Pi*v); % sinus, freq=v
figure;
plot(x,y);

%% filter def
vc_val=1; %1000 Hz
vccont=1/Pi/Ts*tan(Pi*Ts*vc_val);

a=Pi^2*Ts^2*vccont^2;
b=sqrt(2)*Pi*Ts*vccont;

%% filtering
yf=zeros(1,size(y,2));

y1=0; y2=0; y3=0; y4=0; x1=0; x2=0; x3=0; x4=0;

for i=3:size(y,2)
     if (i>1)
         y1=yf(i-1);
         x1=xyi-1);
         if (i>2)
             y2=yf(i-2);
             x2=y(i-2);
         end
     end
    yf(i)=1/(a-b+1)*(-(2*a-2)*y1-(a+b+1)*y2+y(i)-2*x1+x2);
end

figure;
plot(x,yf);

但它并没有给我我期望的结果（左：原始窦，右：切割频率的 1/2 处的过滤结果）：

我检查了我的过滤器是否稳定，我看不出可能是什么问题。你有什么提示吗？

谢谢！

Answer 1

结果可能是正确的。您的正弦波频率为 500，采样频率为 10000。您的归一化频率为 500/10000=0.05，已被过滤。所以你会看到高度衰减的输出。你应该增加你的正弦波频率，看看你是否在过滤后恢复了你的正弦波。

Answer 2

我发现了问题，就是matlab的numden()函数没有给出正确的分子和分母。感谢那些花时间帮助我的人！

以下代码给出了错误的结果：

z=sym('z');
Pi=sym('Pi');
Ts=sym('Ts');
vc=sym('vc');

Fz=1/((Pi^2*Ts^2*vc^2)/(z^2 - 2*z + 1) + (2^(1/2)*Pi*Ts*vc)/(z - 1) + (2*Pi^2*Ts^2*vc^2*z)/(z^2 - 2*z + 1) + (Pi^2*Ts^2*vc^2*z^2)/(z^2 - 2*z + 1) + (2^(1/2)*Pi*Ts*vc*z)/(z - 1) + 1)

[num, den]=numden(Fz);
cn=coeffs(num,z)
cd=coeffs(den,z)

matlab 给出系数：

cn =

[ 1, -2, 1]

cd =

[ Pi^2*Ts^2*vc^2 - 2^(1/2)*Pi*Ts*vc + 1, 2*Pi^2*Ts^2*vc^2 - 2, Pi^ 2*Ts^2*vc^2 + 2^(1/2)*Pi*Ts*vc + 1]

明显是错误的

谢谢！