综述

前情提要

本文使用LMS自适应滤波器算法进行电赛G题的发挥部分的解算。所有算法均在时域中解算,不涉及任何频域!(被FFT干烂的同学不要害怕)
本文默认读者能够完成本题的基础要求部分,以及基础的数字信号处理知识。(交流时候发现有人不知道FIR有限长数字滤波器是怎么滤波的,给整了个大无语)
本人是负责硬件部分的,但是因为学期末数字信号处理的课设项目涉及到了这个算法,所以在此故意卖弄。我有的代码都已经贴在这篇文章里面了,要别的代码我是一点没有的。

本文结构【省流模式】:

  • 2025年电赛G题题目:给出了原题
  • LMS自适应滤波器:给出了算法原理和使用方案,不想看理论,可以直接从改造篇章开始看。
  • MATLAB仿真:结合仿真探讨如何设计参数
  • 可行方案:给出了硬件搭建的结构
  • 滤波器性质判断:解决了原本看起来最简单,但是需要拐个弯的东西
  • 附录:给出了完整的MATLAB仿真代码

2025年电赛G题题目:

任务

设计并制作RC有源滤波电路(简称已知模型电路)和电路模型探究装置(简称探究装置)。基本要求中,探究装置可自动调节本身的输出信号并加给“已知模型电路”,使该电路按要求输出信号;发挥部分中,探究装置可对测评现场提供的“未知模型电路”进行自主学习、建模,并根据该电路输入端的信号推理生成与该电路相同的输出信号。

要求

  1. 基本要求
    (1) 在单独的电路板上搭建电压传递函数为 H ( s ) = 5 1 0 − 8 s 2 + 3 × 1 0 − 4 s + 1 H(s)=\frac{5}{10^{-8}s^2+3\times10^{-4}s+1} H(s)=108s2+3×104s+15的“已知模型电路”。由信号发生器输入频率100Hz~3kHz、峰峰值 1V 的正弦信号,使用示波器测量该电路输出电压幅度,如图1所示。要求该电路输出电压幅度与传递函数 H ( s ) H(s) H(s)表征的输出电压幅度的相对误差绝对值不大于10%。
    (2)探究装置能产生正弦信号,要求频率可设置(步长100Hz),最高频率不小于1MHz,频率相对误差绝对值不大于5%;各频点输出电压峰峰值的最大值不小于3V。
    (3)将探究装置的输出端连接“已知模型电路”的输入端,如图2所示。设置探究装置输出1kHz的正弦信号,并依据基本要求(1)中的H(s)确定输出信号幅度,使得“已知模型电路”输出电压峰峰值为2V。要求该电路输出电压与设定值(2V)的相对误差绝对值不大于5%。
    (4)基本要求(3)中,探究装置可设置并输出100Hz ~ 3kHz频率范围内的正弦信号,依据基本要求(1)中的 H ( s ) H(s) H(s)确定输出信号幅度,使得“已知模型电路”输出电压峰峰值为设定值(范围1 ~ 2V、步长0.1V)。要求该电路输出电压与设定值的相对误差绝对值不大于5%。
  2. 发挥部分
    (1)测评现场提供由RLC元件(各1个)组成的“未知模型电路”。按照图3所示,探究装置连接该电路的输入和输出端口,对该电路进行自主学习、建模(不可借助外部测试设备),2分钟内完成学习建模,显示“未知模型电路”的滤波类型。
    图3 学习、建模
    (2)学习建模完成后,断开探究装置和“未知模型电路”的端口连接,将信号发生器接入探究装置和“未知模型电路”的输入端口(探究装置输入电阻不小于100kΩ),如图4所示。探究装置能根据信号发生器的输出信号推理生成与“未知模型电路”相同的输出信号。要求探究装置和“未知模型电路”的输出信号相比波形无失真、在示波器上能连续同频稳定显示(相位无要求),两者峰峰值相对误差绝对值不大于10%。 信号发生器的输出为频率1kHz ~ 50kHz(步长200Hz)、峰峰值2V的周期信号,类型分别为:正弦波、矩形波(占空比10% ~ 50%、步长5%)和其他周期信号。
    在这里插入图片描述

(3)其他。

LMS自适应滤波器

原理

传统的数字滤波的权系数是固定的,在设计阶段选定,并在滤波器的正常运行中保持不变,然而自适应滤波器是利用前一时刻已获得的滤波器参数等,自动地调节、更新现时刻的滤波器参数,以适应信号和噪声未知的统计特性或随时间变化的统计特性,从而实现最优滤波。自适应滤波器的原理方框图如下所示:
在这里插入图片描述
其中, x ( n ) x(n) x(n) 表示自适应滤波器的输入信号, y ( n ) 表示滤波器的输出信号 y(n) 表示滤波器的输出信号 y(n)表示滤波器的输出信号 d ( n ) d(n) d(n)表示滤波器的参考信号, ϵ ( n ) \epsilon(n) ϵ(n) 表示滤波器的误差信号。
自适应滤波器是由参数可调的数字滤波器 F i l t e r ( ⋅ ) Filter(\cdot) Filter()和自适应滤波算法 C o e ( ⋅ ) Coe(\cdot) Coe()构成。 F i l t e r ( ⋅ ) Filter(\cdot) Filter()的本质是一个 N N N位的有限单位冲击响应数字滤波器(FIR),只是在 n + 1 n+1 n+1时刻的权值 ω i ( n + 1 ) \omega_i(n+1) ωi(n+1)将会由 n n n时刻的权值 ω i ( n ) \omega_i(n) ωi(n) n n n时刻自适应滤波算法 C o e ( n ) Coe(n) Coe(n)得到。该滤波器的信号关系为:
{ y ( n ) = F i l t e r [ x ( n ) ] ϵ ( n ) = d ( n ) − y ( n ) w i ( n + 1 ) = w i ( n ) + C o e [ x ( n ) , d ( n ) , ϵ ( n ) ] \begin{cases} y(n)=Filter[x(n)]\\ \epsilon(n)=d(n)-y(n)\\ w_i(n+1)=w_i(n)+Coe[x(n),d(n),\epsilon(n)] \end{cases} y(n)=Filter[x(n)]ϵ(n)=d(n)y(n)wi(n+1)=wi(n)+Coe[x(n),d(n),ϵ(n)]
由于 F i l t e r ( ⋅ ) Filter(\cdot) Filter()是FIR型滤波器,此处可直接写出滤波器的传递关系:
y ( n ) = F i l t e r ( x ( n ) ) = W ⊤ ( n ) X ( n ) = ∑ i = 0 N − 1 w i ( n ) x ( n − i ) y(n)=Filter(x(n))=\bm{W}^\top(n)\bm{X}(n)=\sum_{i=0}^{N-1}w_i(n)x(n-i) y(n)=Filter(x(n))=W(n)X(n)=i=0N1wi(n)x(ni)
滤波器的目标是在尽可能去除噪声的同时,尽可能的保留原信号。
若将含噪的信号 d ( n ) d(n) d(n),视作希望被提取信号 s ( n ) s(n) s(n)和噪声信号 x ′ ( n ) x'(n) x(n)的线性叠加。根据最小均方误差的思想可以得到第 k k k个点的代价函数为:
C ( k ) = ( Y − S ) ⊤ ( Y − S ) = ∑ i = 0 k ( y ( i ) − s ( i ) ) C(k)=\sqrt{(\bm{Y}-\bm{S})^\top(\bm{Y}-\bm{S})}=\sqrt{\sum_{i=0}^{k}(y(i)-s(i))} C(k)=(YS)(YS) =i=0k(y(i)s(i))
显然,代价函数越小,表明噪声去除得越干净。
由于实际情况中,我们难以完美的得到理想的提取信号 s ( n ) s(n) s(n),因此我们希望通过对输入的噪声信号进行最优估计,也就是滤波器的输出 y ( n ) y(n) y(n)以此来得到提取信号 s ( n ) s(n) s(n)的最优估计 ϵ ( n ) \epsilon(n) ϵ(n)。系统误差 ϵ \epsilon ϵ为:
ϵ = s + x ′ − y \epsilon=s+x'-y ϵ=s+xy
将上式两边平方后得:
ϵ 2 = s 2 + ( x ′ − y ) 2 + 2 s ( x ′ − y ) \epsilon^2=s^2+(x'-y)^2+2s(x'-y) ϵ2=s2+(xy)2+2s(xy)
各取数学期望值,因为 s s s x ′ x' x y y y不相关,因此 E [ s ( x ′ − y ) ] = 0 \mathbb{E}[s(x'-y)]=0 E[s(xy)]=0,所以有:
E [ ϵ 2 ] = E [ s 2 ] + E [ ( x ′ − y ) 2 ] + 2 E [ s ( x ′ − y ) ] = E [ s 2 ] + E [ ( x ′ − y ) 2 ] \mathbb{E}[\epsilon^2]=\mathbb{E}[s^2]+\mathbb{E}[(x'-y)^2]+2\mathbb{E}[s(x'-y)]=\mathbb{E}[s^2]+\mathbb{E}[(x'-y)^2] E[ϵ2]=E[s2]+E[(xy)2]+2E[s(xy)]=E[s2]+E[(xy)2]
当条件滤波参数使得 E [ ϵ 2 ] \mathbb{E}[\epsilon^2] E[ϵ2]最小化时,不希望有效信号能量 E [ s 2 ] \mathbb{E}[s^2] E[s2]受到损失,从而就应当使 E [ ( x ′ − y ) 2 ] \mathbb{E}[(x'-y)^2] E[(xy)2]项最小化,即:
E m i n [ ϵ 2 ] = E [ s 2 ] + E m i n [ ( x ′ − y ) 2 ] \mathbb{E}_{min}[\epsilon^2]=\mathbb{E}[s^2]+\mathbb{E}_{min}[(x'-y)^2] Emin[ϵ2]=E[s2]+Emin[(xy)2]
因此滤波器的输出 y y y应当是噪声信号 x ′ x' x的最小估计。滤波器输出信号 y y y是输入信号 x x x的滤波值。所以,只有输入信号 x x x和噪声信号 x ′ x' x相关,才能有效的使误差的信号数学期望变小。
得误差序列的均方值 Γ \Gamma Γ为:
Γ = E [ ϵ 2 ( n ) ] = E [ ( d ( n ) − y ( n ) ) 2 ] \Gamma=\mathbb{E}[\epsilon^2(n)]=\mathbb{E}[(d(n)-y(n))^2] Γ=E[ϵ2(n)]=E[(d(n)y(n))2]
将此式代入滤波器传递函数有:
Γ = E [ d 2 ( n ) ] + W ⊤ ( n ) R W ( n ) − 2 W ⊤ ( n ) P \Gamma=\mathbb{E}[d^2(n)]+\bm{W}^\top(n)\bm{R}\bm{W}(n)-2\bm{W}^\top(n)\bm{P} Γ=E[d2(n)]+W(n)RW(n)2W(n)P
其中, R \bm{R} R N N N阶自相关方阵,表示输入信号采样值之间的相关性矩阵; P = E [ d ( n ) X ( n ) ] \bm{P}=\mathbb{E}[d(n)\bm{X}(n)] P=E[d(n)X(n)] N × 1 N\times1 N×1互相关矩阵,表示理想信号 d ( n ) d(n) d(n)与输入信号矢量的相关性。
在均方差 Γ \Gamma Γ最小时,最佳权系数 W ∗ = [ w 0 ∗ , w 1 ∗ , ⋯   , w N − 1 ∗ ] \bm{W}^*=[w_0^*,w_1^*,\cdots, w_{N-1}^*] W=[w0,w1,,wN1]应满足梯度向量为零向量:
∇ ( n ) = ∂ Γ ∂ W ( n ) ∣ W ( n ) = W ∗ = 0 \bm{\nabla}(n)=\left.\frac{\partial \Gamma}{\partial \bm{W}(n)}\right|_{\bm{W}(n)=\bm{W}^{*}}=\bm{0} (n)=W(n)Γ W(n)=W=0
即,
R W ∗ − P = 0 \bm{R}\bm{W}^*-\bm{P}=\bm{0} RWP=0
对于这个线性方程组,如果 R \bm{R} R满秩,则有逆矩阵 R − 1 \bm{R}^{-1} R1存在,可以得到权系数最佳值:
W ∗ = R − 1 P \bm{W}^*=\bm{R}^{-1}\bm{P} W=R1P
该方程是 Wiener-Hopf 方程的矩阵形式。
LMS自适应算法是一种实时查找的近似解的实用方法。根据这种方法,下一时刻的权重向量 W ( n + 1 ) W(n+1) W(n+1)等于当前时刻权重向量 W ( n ) W(n) W(n)加上与负梯度成正比的变化,即公式:
W ( n + 1 ) = W ( n ) − μ ∇ ( n ) \bm{W}(n+1)=\bm{W}(n)-\mu\bm{\nabla}(n) W(n+1)=W(n)μ(n)
式中, μ \mu μ表示自适应步长, ∇ ( n ) \bm{\nabla}(n) (n) n n n次迭代的梯度,表示为:
∇ ( n ) = ∂ E [ ϵ 2 ( n ) ] ∂ W ( n ) \bm{\nabla}(n)=\frac{\partial \mathbb{E}[\epsilon^2(n)]}{\partial \bm{W}(n)} (n)=W(n)E[ϵ2(n)]
直接使用瞬时响应替代数学期望,可以得到:
∇ ^ ( n ) = ∂ ϵ 2 ( n ) ∂ W ( n ) = − 2 ϵ ( n ) X ( n ) \bm{\hat{\nabla}}(n)=\frac{\partial \epsilon^2(n)}{\partial \bm{W}(n)}=-2\epsilon(n)\bm{X}(n) ^(n)=W(n)ϵ2(n)=2ϵ(n)X(n)
将瞬时响应下的梯度代入式LMS算法即可得到自适应滤波器的滤波器系数梯度下降算法公式:
W ( n + 1 ) = W ( n ) + 2 μ ϵ ( n ) X ( n ) \bm{W}(n+1)=\bm{W}(n)+2\mu\epsilon(n)\bm{X}(n) W(n+1)=W(n)+2μϵ(n)X(n)

改造

根据上述理论可以知道,自适应滤波器本质是,将与噪声相关的信号通过参数可调的FIR数字滤波器 F i l t e r ( ⋅ ) Filter(\cdot) Filter()拟合含噪信号中的噪声信号部分,此时两信号做减法,便去除了噪声信号,提取得到有用的信号 ϵ ( n ) \epsilon(n) ϵ(n),换句话说FIR滤波器 F i l t e r ( ⋅ ) Filter(\cdot) Filter()的输出 y ( n ) y(n) y(n)不是滤除噪声后的有效信号,而是拟合 d ( n ) d(n) d(n)中的噪声部分。(如果原理部分跳过或没看懂,请务必对照自适应滤波器的原理方框图理解这句话)

我们大胆假设“含噪信号” d ( n ) d(n) d(n)是中有用的信号为0,那么滤波器输出信号 y ( n ) y(n) y(n),就是对参考信号 d ( n ) d(n) d(n)的最优拟合,即 y ( n ) ≈ d ( n ) y(n)\approx d(n) y(n)d(n)。这不就是题目要求我干的探究装置拟合未知模型电路嘛?

于是我们得到了下面的这张示意图:
在这里插入图片描述

  • 绿线和紫线是处理器内部的信号传递,黑色是模拟信号的传递。
  • 蓝色是探究装置整体方框图,包含数据处理模块(如MCU、FPGA、DSP),DA/AD模块两样,是完成提高部分的必要装置。
    = 紫色是自适应学习算法,本质是梯度下降求局部最优的FIR数字滤波器参数。
  • 虚线框中的相移是指因为计算或者输出滞后而导致的相位偏差,这里的相位偏差对于不同的信号通路是不同的,且不可计算得到单独一路的绝对相偏。但是我们可以保证两路的相对相位偏差是一致的,所以,我们建立的LMS自适应滤波器本质上是对于“未知模型电路+未知相位的移相器”的最优估计。

提高部分第一问,摘去示波器,添加学习部分,即可自主学习;第二问,摘去学习部分,通过DAC将 y ( n ) y(n) y(n)输出出来,即为模仿输出。(未知相位的移相器会使得信号相位产生偏差,但题目恰恰说明相位无要求)

滤波器公式

设学习率 μ \mu μ为常数 μ 0 \mu_0 μ0(通常 0 < μ 0 ≤ 0.5 0<\mu_0\leq0.5 0<μ00.5), k k k阶FIR数字滤波器在 n n n时刻的权系数为 w i , n ( i ∈ 0 , 1 , ⋯   , k − 1 ) w_{i,n}(i\in{0,1,\cdots,k-1)} wi,n(i0,1,,k1) n n n时刻的输入信号为 x ( n ) x(n) x(n),未知模型电路对应输出值为 d ( n ) d(n) d(n)
{ y ( n ) = ∑ i = 0 k − 1 w i , k ( n ) x ( n − i ) ϵ ( n ) = d ( n ) − y ( n ) w i , k + 1 = w i , k + 2 μ 0 ϵ ( n ) x ( n − i ) \begin{cases} y(n)=\sum_{i=0}^{k-1}w_{i,k}(n)x(n-i)\\ \epsilon(n)=d(n)-y(n)\\ w_{i,k+1}=w_{i,k}+2\mu_0\epsilon(n)x(n-i) \end{cases} y(n)=i=0k1wi,k(n)x(ni)ϵ(n)=d(n)y(n)wi,k+1=wi,k+2μ0ϵ(n)x(ni)
以上便是学习的全部过程。第二问只需将后两个公式删去,保留第一个公式作为输出即可。

使用注意事项

  1. 因为参数可调的数字滤波器是未知模型电路的最优拟合,所以理论上不管学习还是模仿输出,输入信号均无特殊要求,白噪声信号都行。但在实际中,会受到AD/DA的带宽限制,所以建议还是使用正弦低频信号。(不要自找麻烦)
  2. 因为所有计算、学习均在时域内拟合完成,因此第一问的学习部分重要的是采样点个数,而不是采样频率。只要采样频率不是零理论上都能学收敛,但是两分钟的限时肯定完蛋。
  3. 阶数太大会使得收敛效率很低,两分钟可能跑不完;阶数太小会使得收敛结果并能很好的拟合未知模型电路。这里推荐阶数不低于4阶。
  4. 你会发现,第二问其实是阉割的第一问。所以这个方案的风险性巨大,要么全成功,要么全G!

MATLAB仿真

核心代码

以下是LMS自适应滤波器学习算法的MATLAB代码。全部MATLAB代码在附录中给出。
我只试过32阶和10阶;学习率0.01和0.001的仿真,对应50000-100000点就能学出很好的模样(10kHz采样率,5-10s采完)。

% 输入信号x(n)是代码中的数组u
% 未知模型电路对应输出值d(n)是代码中的数组y_lagged

% 定义滤波器参数
filter_order = 10;  % 滤波器长度
w = ones(filter_order, 1); % 初始权重(归一化)
mu = 0.01;        % 学习率(减小以避免发散)

% 初始化输出
p = zeros(size(u)); % 滤波后信号
q = zeros(size(u)); % 误差信号

% 自适应滤波处理
for k = filter_order:length(u)-1
    % 获取当前输入窗口
    current_window = u(k-filter_order+1:k);
    
    % 计算滤波输出
    p(k) = w' * current_window;
    
    % 计算误差(注意滞后对齐)
    if k > lag_points
        q(k) = p(k) - y_lagged(k);
    else
        q(k) = 0; % 初始阶段没有滞后信号可用
    end
    
    % LMS权重更新
    w = w - 2*mu*q(k)*current_window;
end

仿真现象及分析

仿真使用的信号为100Hz到3kHz正弦波信号,随机频率,每200点更换频率值,直至采样结束。

阶数10,学习率0.01,采样率10kHz,采样时长0.1s

在这里插入图片描述
为了方便观察图像讨论,此处调低了采样学习时长,数字滤波器并不能良好的复刻未知模型输出,所以最下方的模仿输出显然是不正确的。

上图主要看绿色曲线,绿色曲线是误差值 ϵ ( n ) = y ( n ) − d ( n ) \epsilon(n)=y(n)-d(n) ϵ(n)=y(n)d(n)

  • 对于单个频率的正弦信号,200点足够自适应滤波器学习收敛。所以扩展单频率的采样点数是没有意义的。
  • 对于多个频率的交界处,这个地方是最容易出现问题的地方。此处涉及频率的交接,频谱分量丰富,是数字滤波器最主要的学习之处。

因此,对于采样时长和学习效果的优化,最理想的方案是在单频率信号稳定跟随时立刻切换频率值,以此加快学习进程。

阶数10,学习率0.001,采样率10kHz,采样时长30s

在这里插入图片描述
对于10阶滤波器,我按标题给出的参数跑了100遍仿真,有31遍不满足峰峰值相对误差绝对值不大于10%的要求,接近三分之一的概率。
以上是其中一遍仿真的输出信号,可以看到这轮输出由于相位是对上了所以很容易看到除了幅值不足以外,相位,信号模样基本上都是一模一样的。这也说明了学习的成功性。

阶数32,学习率0.01,采样率10kHz,采样时长10s

在这里插入图片描述
32阶FIR已经可以说是完美复刻了,从4秒起误差已经基本上没有了波动,稳定在0附近。

在100次学习仿真后全部都符合题目的幅值要求。这说明,阶数不足的时候,再怎么学都没有用。必须提升阶数,才能更好的复刻。

可行方案

比赛采用的方案

我们比赛用正点原子的FPGA作为主控制,负责信号生成和学习未知模型电路。只将低通,高通,带通,带阻判断结果发送给MCU用作屏幕输出。主框图如下所示:
在这里插入图片描述
但是,我们发现FPGA写LMS算法十分困难,所以不推荐用FPGA来做这道题。
但是这里可以贴一张比赛时候的调试图:
在这里插入图片描述
低通滤波器的效果已经学出来了。(这似乎是我唯一能找到的图)

似乎更好的方案

在FPGA队友库库写代码的时候,我用STM32H723试了一下,感觉用MCU代替FPGA会更加方便(使用专门的并口AD/DA模块能够做到±5V输入输出,这个应当比32自带的AD/DA做题方便)。

对应LMS自适应滤波器的C语言代码非常非常非常非常非常非常非常非常非常非常简单,核心代码如下:

#ifndef __LMS_H__
#define __LMS_H__

#define FILTER_LEN 10 //滤波器长度

typedef struct{
	float w[FILTER_LEN];
	float adcin[FILTER_LEN];
	float error;
	float mu; 
}LMS_struct;

void  LMS_init       (LMS_struct& LMS,float mu);
float LMS_study_cul  (LMS_struct& LMS,float dds_value,float filter_value);
float LMS_output_cul (LMS_struct& LMS,float dds_value);
#endif
#include "LMS.h"
//初始化,mu为学习率
void LMS_init(LMS_struct& LMS,float mu){
	for(int i=0;i<FILTER_LEN;i++){
		LMS->w[i]=1;
		LMS->adcin[i]=0;
	}
	LMS->error=0;
	LMS->mu=mu;
}

//学习LCR,dds_value是x(n)的采样值,filter_value是d(n)的采样值
//函数输出数字滤波器输出y(n)
float LMS_study_cul(LMS_struct& LMS,float dds_value,float filter_value){
	float output=0;
	for(int i=FILTER_LEN-1;i>0;i--){LMS->adcin[i]=LMS->adcin[i-1];}
	LMS->adcin[0]=dds_value;
	for(int i=0;i<FILTER_LEN;i++){output+=LMS->adcin[i]*LMS->w[i];}
	LMS->error=output-filter_value;
	for(int i=0;i<FILTER_LEN;i++){LMS->w[i]+=2*mu*LMS->error*LMS->adcin[i];}
	return output;
}

//复刻LCR滤波输出,dds_value是x(n)的采样值
//函数输出数字滤波器输出y(n)
float LMS_output_cul (LMS_struct& LMS,float dds_value){
	float output=0;
	for(int i=FILTER_LEN-1;i>0;i--){LMS->adcin[i]=LMS->adcin[i-1];}
	LMS->adcin[0]=dds_value;
	for(int i=0;i<FILTER_LEN;i++){output+=LMS->adcin[i]*LMS->w[i];}
	return output;
}

滤波器性质判断

关于这一点,我们比赛是没做的,因为FPGA实在是太难写了(FPGA队友原话)。如果用MCU做那就是信手拈来了。

事先准备好多个频率的200点采样值,然后用学习好的FIR数字滤波器去算,取出每个频率信号的最大值作为幅值直接比较就可以了得到大致的伯德图,相当于是扫频,但这个扫频完全是数学卷积运算,应当有手就行。

附录

完整的MATLAB仿真代码复制如下,如果有需要的同学可以直接使用

% Apercent是
function Apercent=filttttter()
    % 二阶传递函数对正弦信号的响应分析
    clear all;
    close all;
    clc;
    
    % 创建传递函数
    num = 1;
    den = [2.2e-11,1e-5,1];
    sys = tf(num, den);
    
    % 基本参数设置
    fs = 10000;              % 采样频率(Hz)
    total_time = 10;         % 总时长(s)
    num_points = fs * total_time; % 总点数
    t = (0:num_points-1)/fs; % 时间向量
    
    % 分段设置(每50个点为一个段)
    segment_length = 200;
    num_segments = floor(num_points/segment_length);
    
    % 初始化频率数组(每个段一个随机频率)
    min_freq = 100;          % 最小频率(Hz)
    max_freq = 3000;         % 最大频率(Hz)
    frequencies = min_freq + (max_freq-min_freq)*rand(1, num_segments);
    
    % 生成变频正弦信号
    signal = zeros(1, num_points);
    for seg = 1:num_segments
        start_idx = (seg-1)*segment_length + 1;
        end_idx = min(seg*segment_length, num_points);
        
        % 当前段的频率
        current_freq = frequencies(seg);
        
        % 生成当前段的正弦信号
        signal(start_idx:end_idx) = sin(2*pi*current_freq*t(start_idx:end_idx));
    end
    
    u=signal';
    
    % 计算系统响应
    [y, t] = lsim(sys, u, t);
    
    % 滞后输出信号
    lag_points = 5; % 滞后点数
    y_lagged = [zeros(lag_points, 1); y(1:end-lag_points)];
    
    % 定义滤波器参数
    filter_order = 32;  % 滤波器长度
    w = ones(filter_order, 1); % 初始权重(归一化)
    mu = 0.01;        % 学习率(减小以避免发散)
    
    % 初始化输出
    p = zeros(size(u)); % 滤波后信号
    q = zeros(size(u)); % 误差信号
    
    % 自适应滤波处理
    for k = filter_order:length(u)-1
        % 获取当前输入窗口
        current_window = u(k-filter_order+1:k);
        
        % 计算滤波输出
        p(k) = w' * current_window;
        
        % 计算误差(注意滞后对齐)
        if k > lag_points
            q(k) = p(k) - y_lagged(k);
        else
            q(k) = 0; % 初始阶段没有滞后信号可用
        end
        
        % LMS权重更新
        w = w - 2*mu*q(k)*current_window;
    end
    
    subplot(4,1,1);
    plot(t, u, 'k', 'LineWidth', 1);
    title('输入信号(带噪声)');
    ylabel('幅值');
    grid on;
    
    subplot(4,1,2);
    plot(t, y, 'b', 'LineWidth', 1.5); 
    hold on;
    plot(t, p, 'r', 'LineWidth', 1.5);
    title('系统响应 vs 自适应滤波输出');
    legend('系统响应', '自适应滤波');
    ylabel('幅值');
    grid on;
    
    subplot(4,1,3);
    plot(t, q, 'g', 'LineWidth', 1);
    title('误差信号');
    xlabel('时间 (s)');
    ylabel('误差');
    grid on;
    
    % 输入正弦信号参数
    A = 1;           % 幅值 (固定为1)
    qq = 290*randn(1)+1000;           % 频率 (rad/s)
    phi = pi/4;      % 相位 (rad)
    
    % 仿真时间设置
    t = linspace(0, 10, 1000)'; % 时间向量(列向量)
    
    % 生成输入信号
    pwm_signal = sin(2*pi*qq*t + phi);
    
    % % 参数设置
    % fs = 100000;       % 采样频率 (Hz),需至少是信号频率的10倍(Nyquist定理)
    % f = 10000;         % 方波频率 (10 kHz)
    % duration = 0.01;   % 信号持续时间 (秒)
    % t = 0:1/fs:duration; % 时间向量
    % 
    % % 生成方波(占空比50%)
    % duty_cycle = 50;   % 占空比(50%为对称方波)
    % pwm_signal = (0.5 * (1 + sign(sin(2 * pi * f * t))))';
    
    % 数字滤波处理
    filttter = zeros(size(pwm_signal)); % 滤波后信号
    for k = filter_order:length(pwm_signal)-1
        % 获取当前输入窗口
        current_window = pwm_signal(k-filter_order+1:k);
        
        % 计算滤波输出
        filttter(k) = w' * current_window;
    end
    
    % 计算系统响应
    [system_response, t_response] = lsim(sys, pwm_signal, t);
    
    subplot(4,1,4);
    plot(t, system_response, 'b', 'LineWidth', 0.5);
    hold on;
    plot(t, filttter, 'r', 'LineWidth', 0.5);
    title('系统响应 vs 模拟输出');
    xlabel('时间 (s)');
    ylabel('幅值');
    
    Apercent=(1-max(filttter)/max(system_response))*100;
end
Logo

智能硬件社区聚焦AI智能硬件技术生态,汇聚嵌入式AI、物联网硬件开发者,打造交流分享平台,同步全国赛事资讯、开展 OPC 核心人才招募,助力技术落地与开发者成长。

更多推荐