如何用MATLAB的FDA/FilterDesigner设计低通滤波器保留±40π信号？

我正在为信号D(ω)设计滤波器，目标是移除4960π和5040π处的干扰信号，仅保留±40π处的有效信号。考虑到Chebyshev滤波器过渡带更窄，选用MATLAB内置的FilterDesigner工具进行设计（Butterworth滤波器也可），但运行MATLAB脚本时出现错误，且无图像输出。

相关截图 #

完整MATLAB代码 #

%create symbolic functions x, a, b, c, d, e, f_c with independent variable t
syms x(t) a(t) h(t) b(t) c(t) d(t) e f_c(t) f_c1(t) f_c2(t) t tau

%create symbolic functions A, B, C, D, and E with independent variable w
syms A(w) B(w) C(w) D(w) E(w) w

x(t) = cos(100*pi*t);
a(t) = x(0.4*t);
h(t) = dirac(t-0.02);
b(t) = int(a(tau)*h(t-tau), 'tau', -inf, inf);

f_c(t) = 10*cos(2500*pi*t);
f_c1(t) = f_c(t);
f_c2(t) = f_c(t);

c(t) = b(t)*f_c1(t);
d(t) = c(t)*f_c2(t);

figure

subplot (3,1,1)
fplot(a(t))
xlim([-0.05 0.05]),ylim([-1.5 1.5])
title ('Time domain of signal a(t)')
xlabel('Time, t')
ylabel('Amplitude, a(t)')
grid on

A(w) = fourier(a(t), w);
w = -60*pi:0.1*pi:60*pi;
subsA = A(w);
% Replace Inf value with suitable value
idx = subsA == Inf;
subsA(idx) = pi; % choose suitable value from the expression of fourier transform.
subplot (3,1,2)
plot(w,real(subsA));
ylim([0 4])
title ('Frequency domain of signal A(\omega)')
ylabel("\Re(A(\omega))");
xlabel("\omega");
xticks(-40*pi:20*pi:40*pi);
xticklabels({'-40\pi', '-20\pi', '0', '20\pi', '40\pi'})

subplot (3,1,3)
plot(w, angle(A(w)))
title ('Phase spectrum of signal A(\omega)')
ylabel("\angle(A(\omega))");
xlabel("\omega");


figure

subplot(3,1,1)
fplot(b(t))
xlim([-0.05 0.05]),ylim([-1.5 1.5])
title ('Time domain of signal b(t)')
xlabel('Time, t')
ylabel('Amplitude, b(t)')
grid on

syms B w
B(w) = fourier(b(t), w);
w = -60*pi:0.1*pi:60*pi;
subsB = B(w);
% Replace Inf value with suitable value
idx = abs(subsB) == Inf;
subsB(idx) = pi; % choose suitable value from the expression of fourier transform.
subplot (3,1,2)
plot(w,real(subsB));
ylim([0 4])
title ('Frequency domain of signal B(\omega)')
ylabel("\Re(B(\omega))");
xlabel("\omega");
xticks(-40*pi:20*pi:40*pi);
xticklabels({'-40\pi', '-20\pi', '0', '20\pi', '40\pi'})

subplot (3,1,3)
plot(w, angle(B(w)))
ylim([-4 4])
title ('Phase spectrum of signal B(\omega)')
ylabel("\angle(B(\omega))");
xlabel("\omega");

figure

subplot(3,1,1)
fplot(c(t))
xlim([-0.05 0.05]),ylim([-12 12])
title ('Time domain of signal c(t)')
xlabel('Time, t')
ylabel('Amplitude, c(t)')
grid on

syms C w
C(w) = fourier(c(t), w);
%simplify(C(w))
w = -2800*pi:20*pi:2800*pi;
subsC = C(w);
subsC
% Replace Inf value with suitable value
idx = abs(subsC) == Inf;
idx
subsC(idx) = pi; % choose suitable value from the expression of fourier transform.
subplot (3,1,2)
plot(w,real(5*subsC));
%testing = subsC
ylim([0 20])
title ('Frequency domain of signal C(\omega)')
ylabel("\Re(C(\omega))");
xlabel("\omega");


subplot (3,1,3)
plot(w, angle(C(w)))
ylim([-4 4])
title ('Phase spectrum of signal C(\omega)')
ylabel("\angle(C(\omega))");
xlabel("\omega");


figure 

subplot(3,1,1)
fplot(d(t))
xlim([-0.1 0.1]),ylim([-110 110])
title ('Time domain of signal d(t)')
xlabel('Time, t')
ylabel('Amplitude, d(t)')
grid on

syms D w
D(w) = fourier(d(t), w);
w = -5050*pi:10*pi:5050*pi;
%simplify(D(w))
subsD = D(w);
%D(w)
% Replace Inf value with suitable value
idx = abs(subsD) == Inf;
subsD(idx) = pi; % choose suitable value from the expression of fourier transform.
subplot (3,1,2)
plot(w,real(25*subsD));
ylim([0 160])
title ('Frequency domain of signal D(\omega)')
ylabel("\Re(D(\omega))");
xlabel("\omega");

subplot (3,1,3)
plot(w, angle(D(w)))
ylim([-4 4])
title ('Phase spectrum of signal D(\omega)')
ylabel("\angle(D(\omega))");
xlabel("\omega");

figure

LPF_fil = lpf;
E = filter(LPF_fil, D);
subplot(3,1,1)
plot(w, E)

% fc = 150000;
% fs = 1000000;
% [e_mag,e_ang] = butter(6,fc/(fs/2));
% freqz(e_mag,e_ang)
% e = filter(e_mag,e_ang,d);
% eF=fftshift(fft(e));
% eFm=abs(eF);
% eFa=angle(eF);
% eFa(abs(eF)<T)=0;

核心问题定位 #

• 符号与数值变量冲突：多次重复定义syms B w、syms C w、syms D w，覆盖了之前的数值数组w，导致后续计算逻辑混乱。
• 滤波器调用错误：LPF_fil = lpf;中lpf未定义，且filter函数要求输入时域信号，而非频域符号表达式D。
• 符号傅里叶变换数值化缺陷：直接绘制符号傅里叶变换结果，存在数值不稳定（Inf值处理粗糙）、运算效率低的问题。

修复方案 #

1. 清理重复变量定义 #

删除所有重复的syms B w、syms C w、syms D w语句，确保w作为数值数组时不被重新定义为符号变量。

2. 修正滤波器调用逻辑 #

• 通过FilterDesigner导出滤波器系数（如[b,a]格式），替换未定义的lpf；
• 先将符号信号d(t)转换为数值时域序列，再用filter函数处理：

  t = -0.1:1e-4:0.1;
  d_num = subs(d(t), t, t);
  

3. 替换符号傅里叶变换为数值FFT #

符号运算适合理论推导，数值计算建议使用FFT提升效率与稳定性：

N = length(t);
f = (-N/2:N/2-1)*(1/(max(t)-min(t)));
D_fft = fftshift(fft(d_num));
D_mag = abs(D_fft);

4. 完整滤波处理代码片段 #

% 导入Chebyshev低通滤波器系数（从FilterDesigner导出）
[b,a] = cheby1(6, 0.5, 40*pi/(5050*pi)); % 归一化截止频率

% 生成时域信号
t = -0.1:1e-4:0.1;
d_num = subs(d(t), t, t);

% 滤波处理
e_num = filter(b,a, d_num);

% 绘制滤波后结果
E_fft = fftshift(fft(e_num));
E_mag = abs(E_fft);

figure
subplot(2,1,1)
plot(t, e_num)
title('滤波后信号时域波形')
xlabel('时间 t')
ylabel('幅度')

subplot(2,1,2)
plot(f, E_mag)
xlim([-50*pi 50*pi])
title('滤波后信号频谱')
xlabel('\omega')
ylabel('幅度')

内容的提问来源于stack exchange，提问作者Bryson Fernades