ARTICLE DETAIL

资讯详情

深耕网站建设与运营推广的一线实战洞察。

MATLAB FFT滤波实战:从频谱分析到Simulink波形去噪

MATLAB FFT滤波实战:从频谱分析到Simulink波形去噪 1. FFT滤波的整体思路与适用场景做Simulink仿真的人应该都遇到过这种尴尬模型跑出来的波形明明该是一条光滑的正弦曲线示波器里看却是一堆毛刺甚至分不清是信号还是噪声。想滤掉高频干扰又不确定该选巴特沃斯还是切比雪夫滤波器调参数调到怀疑人生。这时候直接用FFT滤波反而是一条最直观、最容易控制的路。所谓FFT滤波本质上就是先把时域波形通过快速傅里叶变换Fast Fourier Transform搬到频率域看看信号的能量集中在哪些频率上然后把不需要的频率成分直接清零或者衰减最后用逆FFTIFFT把结果搬回时域。整个过程不需要设计任何模拟滤波器原型也不需要处理零极点你面对的就是一张频率幅值图想留什么频率、想砍什么频率一目了然。这种频域掩码Spectral Mask式的处理方式特别适合那些信号和噪声在频谱上分得比较开的场景比如工频干扰、机械振动的高频噪声、采集电路里的白噪声毛刺。这个方案解决的核心问题有三个。一是当你只有一段示波器截图或者一堆.mat数据没有任何滤波器设计指标通带、阻带、纹波时可以直接“看着频谱来滤波”。二是当你需要快速试算不同截止频率对波形的真实影响时FFT滤波只需要改一行代码中的频率上下限重算一次FFT和IFFT比重新设计一版滤波器快得多。三是当你需要对整段数据进行批处理逐个查看哪些频率分量值得保留时FFT频谱图就是最好的可视化依据。适合来参考这条路的主要是这几类人在Simulink里做控制算法仿真、需要对输出波形做后处理的工程师做信号处理实验、手里有一堆.mat或者示波器导出波形文件的研究生以及刚学MATLAB不久、想知道FFT除了画频谱图还能怎么用的新手。这篇文章会从数据怎么进MATLAB开始一直讲到怎么把FFT滤波脚本集成到Simulink模型里中间附带我踩过坑的细节尽量让你看完就能直接上手。2. 数据来源示波器波形与.mat文件的读取在做FFT滤波之前第一步其实是搞清楚数据从哪来、在MATLAB里长什么样。很多人卡住不是死在滤波算法上而是死在数据都读不进来或者读进来了变量结构搞不懂。下面这两种来源最典型。2.1 从Simulink示波器导出波形数据Simulink里的Scope模块是大家最常用的观察窗口但Scope显示出来的只是像素点真正要处理数据你得想办法把波形数据从Scope里导出来。有三种方式我常用按推荐程度排第一种用“Scope模块的日志数据”。在旧版Simulink里Scope会默认把输入信号记录到工作区一个名为ScopeData的结构体里类似ScopeData.time和ScopeData.signals.values这样的字段结构。到了新版R2016b之后新版Scope在配置界面里有一个“Logging”选项卡可以设置“Log data to workspace”一般会生成一个logsout的Simulink.SimulationData.Dataset对象。你可以用logsout.getElement(1).Values.Data和logsout.getElement(1).Values.Time把数据拿出来。第二种直接在模型里放置一个To Workspace模块把要观察的信号连到这个模块上命名比如yout设置输出数组名仿真完之后工作区就会多一个yout变量。如果设置的是“Timeseries”格式那yout.Data就是数据yout.Time就是时间轴。这种方法最可靠因为不依赖Scope的内部机制而且你可以在模型里随意选择要导出的信号。第三种用sim()函数在脚本里跑仿真并指定输出。比如out sim(myModel, StopTime, 10)然后把模型里的输出模块作为返回值再从out里取出信号。这种方式适合批处理一次可以跑几十组不同参数然后自动滤波分析。无论用哪种方式导出后你先用whos和size检查一下变量大小别急着滤波。我试过很多次导出的数据是一列数组但时间轴可能不是等间隔的——特别是变步长仿真后Scope显示的时间点是不均匀的。如果时间轴不等距做FFT之前必须先插值重采样到等间隔否则得到的频谱会有一堆假频谱峰。2.2 外部.mat数据文件的加载与结构解析外部.mat文件就更多样了可能是别的同事发的可能是示波器软件导出的也可能是老版本MATLAB存的。加载之前先用whos(-file, 文件名.mat)看看里面都有什么变量不要一上来就load否则工作区被塞满同名变量覆盖了你的原始数据都不知道。比如文件里存的是一个结构体data里面可能还有字段time和x。典型操作clear; load(wave.mat); whos(-file, wave.mat); % 假设看到 data 是 struct t data.time; x data.x;如果.mat里存的是两个矩阵你想让它们相减再取绝对值那就要先确认尺寸一致。我遇到过很多次矩阵一个是行向量一个是列向量减出来直接报错用size看一眼必要时转置A data.matrixA; B data.matrixB; if size(A) ~ size(B) % 注意size返回向量比较要用 isequal if isequal(size(A), size(B)) % 同构 else % 尝试转置 if isequal(size(A), size(B.)) B B.; else error(矩阵尺寸不匹配); end end end result abs(A - B);另外如果.mat文件里有多个变量你只想要其中一个最好用load(filename, varname)指定加载避免把无关的大数组也读进内存。这一点在数据几十MB以上的时候尤其重要。数据进来了下一步就要开始真正的FFT滤波。但别急先把FFT滤波的基石——参数和预处理——弄明白否则后面全是在瞎猜。3. 核心参数与数据预处理从时间域走向频率域之前很多人直接把原始信号丢进fft()函数画出来的频谱乱七八糟然后抱怨MATLAB的FFT不好用。其实问题不在FFT在于你没有处理好采样的几个常识问题。3.1 采样率、数据长度与频率分辨率之间的硬关系离散信号做傅里叶变换逃脱不了一个铁三角实际采样率Fs、数据点数N、频率分辨率df。它们满足一个极其简单的公式[ df \frac{Fs}{N} ]也就是说你采了1000个点采样率是1000 Hz那频率分辨率就是1 Hz。频率轴第n条谱线对应的实际频率是(n \times df)注意频率范围是0到Fs单边看是0到Fs/2奈奎斯特频率。很多人在这一步犯迷糊。比如Simulink模型用变步长求解器Output的“Save format”选了“Timeseries”但采样率根本不是固定值。如果你直接用fft(x)那FFT默认你的数据间隔是1秒采样率就是1 Hz对应出的频率轴全错。所以在滤波之前你务必确认有没有一个真实的采样周期Ts或者自己指定一个统一的采样时刻。如果原始时间轴不是均匀的我建议先做一次interp1重采样把数据插值到一个统一的Fs上。例如t_uniform t(end) * linspace(0, 1, N_target).; % 生成均匀时间轴 x_uniform interp1(t, x, t_uniform, linear); Fs 1 / (t_uniform(2) - t_uniform(1));这里有个具体的经验重采样点数N_target不要随便拍脑袋它直接影响后面FFT的计算速度和频率分辨率。想让分辨率细一点就增加N_target但分辨率一旦超过你实际采样信息的信息量插值产生的高频细节全是假的滤波时反而会引入奇怪的成分。通常N_target取原始数据点数同等数量级即可比如原始1万点重采样1万点或2万点都行。如果数据本身就是等间隔的那就简单了直接用Fs 1/Ts。这个Ts在你Simulink模型的固定步长里很容易找到在外部采集设备上就是示波器或采集卡的采样间隔。3.2 去直流、去趋势、加窗让频谱干净一半原始信号往往带有直流偏置和趋势项。直流偏置会让频谱0 Hz处有一个大尖峰在进行频域滤波时如果你不是故意保留直流那这个尖峰可能会泄漏到邻近低频淹没真实的低频信号。趋势项就是整个波形缓慢漂移的斜坡对应频谱上很低频的分量同样干扰判断。处理方式很简单先减去均值x x - mean(x);如果你的数据有明显线性漂移还可以先用detrend去掉线性趋势x detrend(x, linear);做完去直流加窗是另外一个容易忽略的细节。直接对一整段数据做FFT相当于在时域上加了一个矩形窗矩形窗的频谱旁瓣很高会导致频率成分向旁边泄漏产生“裙边”效应。尤其在滤除某个窄带强干扰时旁瓣泄漏会让附近频率的幅值也畸变。常规做法是加一个汉宁窗或汉明窗让数据两端渐变为零减少泄漏。不过要注意加窗会改变信号的幅值特别是在恢复时域的绝对值时需要进行幅值修正。窗函数W的等效噪声带宽或相干增益需要校正。通常你可以用w hann(N)然后xw x(:) .* w(:)做FFT滤波后再除以窗均值也就是mean(w)来补偿幅值。很多教程都不提这个导致滤波后波形整体幅度变小还以为滤波有问题。我的建议是先不加窗看一次频谱如果谱线很干净不需要滤除窄带强干扰可以不额外加窗如果需要精确滤除某个窄峰再加上汉宁窗并做幅值补偿。不要盲目加窗。3.3 频率轴的正确生成方式很多初学者喜欢用MATLAB文档自带示例里的freq轴但那个是以采样点数为基数的不方便按实际频率操作。我每次都用这套模板N length(x_uniform); Y fft(x_uniform); f (0:N-1) * (Fs/N); % 双边频率轴然后画幅值谱时通常只画单边f_single f(1:floor(N/2)1); amp abs(Y(1:floor(N/2)1)) / N * 2; % 单边幅值修正直流分量不乘2但这里一般忽略特殊情况 plot(f_single, amp);记住fft输出的第1个点是0 Hz即直流第2个点是Fs/N以此类推。滤波时如果直接对Y做赋值必须确保对称性这一点后面第4节详细说。4. FFT滤波的完整实现从频域掩码到IFFT恢复现在进入正题。假设你有一份干净的等间隔时域数据x采样率Fs接下来我用一段完整代码展示FFT滤波的流程并把里面的坑挑出来讲。4.1 基础版FFT截止滤波场景信号是50 Hz的正弦叠加了5倍幅度的200 Hz噪声想把200 Hz以上滤掉保留50 Hz和附近的低频。代码如下function [x_filt, f, amp] fft_filter_basic(x, Fs, f_low, f_high) % 简单FFT带通滤波保留f_low到f_high之间频率其余置零 N length(x); Y fft(x); % 复数频谱 f (0:N-1) * (Fs/N); % 频率轴 mask ones(N, 1); % 找频率范围对应的索引 idx f f_low | f f_high; mask(idx) 0; % 重要确保频谱共轭对称 Y_filtered Y .* mask; % 如果破坏了共轭对称IFFT会有虚部一般取real x_filt real(ifft(Y_filtered)); amp abs(Y_filtered) / N * 2; end这段代码看着简单但有两个隐患。第一idx会把N个点中所有低于f_low和高于f_high的都置零这确实包括了负频率部分。由于Y是共轭对称的你只置零了正负频率对应的幅值理论上没问题。但如果你不小心只处理了高频部分而没保留对称性——比如自己生成了一个频率索引向量只把正频率部分置零负频率没动——做完IFFT后信号会变成复数产生严重的波形失真。所以请务必记住一句话频域滤波时操作必须同时对正负频率进行或者用完整的频率轴来处理。第二硬截止会在频率边界处产生突变等效于时域乘了一个很长的Sinc函数对应的时域波形会振铃吉布斯效应这在阶跃突变位置尤其明显。如果信号本身是连续光滑的正弦叠加振铃不明显但如果你的信号里有矩形脉冲或阶跃硬截止会在阶跃前后产生明显过冲。这个放在后面第5节专门讲。4.2 带参优化的滤波函数自动识别峰值、保留有效带宽面对实际数据光设一个上下限是不够的。比方说你想滤除某个固定频率的窄带干扰比如50 Hz工频但周围信号能量密集直接整段置零会把附近有用的边频也削掉。这时候更合理的方法是构造一个陷波器或者用高斯型过渡带衰减。我常用的优化版本是允许你指定“保留频率”的峰值点和带宽然后对保留区乘1、禁区乘0、过渡区乘平滑曲线function x_filt fft_filter_smooth(x, Fs, keeps, width) % keeps: 二元矩阵或逻辑向量指示需要保留的频率范围 % width: 过渡带宽度Hz N length(x); Y fft(x); f (0:N-1) * (Fs/N) - Fs/2; % 移动到对称频率轴方便操作 Y_shift fftshift(Y); % 把零频放到中心 % 用一系列高斯过渡带构造掩码 mask ones(N,1); % 先把要保留的中心频率和半径定义出来 for k 1:size(keeps,1) fc keeps(k,1); r keeps(k,2); % 在这个频段内设为1外推高斯衰减 % 简化做法在边界用余弦过渡 end % 伪代码省略实际可用逻辑索引 卷积平滑 end实际做的时候我不会手工写复杂高斯掩码而是先用findpeaks在幅值谱上找到显著峰值再根据峰值周围能量设定保留区间。这样做的好处是避免漏掉那些不起眼的但实际是信号的分量也能避免把噪声峰值误认为信号。大体步骤计算幅值谱amp和频率轴f。用findpeaks(amp, f, MinPeakHeight, ..., MinPeakDistance, ...)找出候选峰值。根据峰值处的半功率带宽-3 dB宽度确定保留区间。除了这些区间外的频率全部清零区间边界用余弦过渡带平滑。这套操作我测过多次对含混叠的采集数据特别有效。4.3 频率变换后注意相位保持很多只关注幅值的教程往往忽略FFT滤波后相位会发生变化。用硬截止滤波保留区相位不变但过了过渡带的频率相位会被置零其实是该频率被删除了时域原有分量没了。这在我们只看幅值时没多大影响但如果信号后续要跟另一个信号做相关分析或求相角那么滤波和被滤之间的相位一致性必须谨慎。IFFT恢复出来的信号跟原始信号在保留频段内相位一致因为只是乘了一个实数掩码没有引入相位偏移。这一点放心。如果使用了平滑过渡带你的掩码是实数向量乘到频谱上也不改变相位。但如果你用了复数掩码比如想改变相频响应那就要小心了。常规FFT滤波我是建议只用实数掩码也就是说每个频率分量要么保留原幅值要么为零/衰减不要动它的相位这样最安全。5. 把FFT滤波接入Simulink模型两种主流方式滤波脚本在MATLAB里跑通了下一步是想办法让它跟Simulink模型配合。这里有两种使用层次我分别说明。5.1 离线后处理导出数据、脚本滤波、回填分析离线方式是最稳妥的。Simulink模型跑完之后把示波器数据导出成yout或者ScopeData然后在MATLAB脚本里执行滤波、画图最后再决定是否把滤波后的信号用在后续分析中。这种方式适合做设计验证、参数扫描、结果报告因为它不干扰仿真本身可以对同一份数据反复调整滤波器参数。举个实际例子。我在调试一个电机控制模型时转速反馈信号里叠加了换相噪声频率大约在300 Hz左右而控制带宽只有50 Hz。模型跑完后我从To Workspace里拿到yout其中yout.signals(1).values是转速yout.time是时间轴。然后t yout.time; x yout.signals(1).values; Fs 1 / median(diff(t)); % 等距时用median更稳 x_filt fft_filter_smooth(x, Fs, [0 60], 10); % 保留0-60Hz过渡带10Hz figure; subplot(2,1,1); plot(t, x); title(原始); subplot(2,1,2); plot(t, x_filt); title(FFT滤波后);这样对比一眼就知道噪声被压下去多少控制周期的相位延迟是否存在。离线处理的好处是随意试错不会导致仿真结果不可重复。如果你想把滤波后的信号送进模型继续仿真那可以把这个滤波结果写回工作区变量然后在下一次仿真时用一个From Workspace模块或Signal Editor导入。注意导入的信号时间轴要和原模型仿真时间匹配否则插值会出问题。5.2 在线集成使用MATLAB Function模块实现实时FFT滤波另一种是直接把FFT滤波搬进Simulink模型里。在模型里放一个MATLAB Function模块双击编辑函数把滤波逻辑写进去输入待滤波信号输出滤波结果。这样每次仿真步长内它都会执行FFT滤波。这里有一个残酷的工程现实FFT是块处理算法天然需要成块的数据而Simulink的离散仿真是一个点一个点算的。如果你每次采样点都单独做一次FFT那个点数通常小得可怜比如1个点根本做不了频域分析。所以在在线方式里你必须维护一个数据缓冲区类似一个滑动窗口。窗口长度取定值N比如1024点。每来一个新的采样点缓冲区滚动更新一次对缓冲区做FFT滤波然后输出最后一个点或者整个窗口。这样有延迟但实时性足够用于离线分析、监控或者预处理。这里推荐两种实现结构第一种使用Interpreted MATLAB Function或MATLAB Function模块加持久变量function y fftfilt_inline(u) persistent buf; if isempty(buf) || true N 1024; buf zeros(N,1); end buf [buf(2:end); u]; % 滑动更新 y fft_filter_basic(buf, Fs, f_low, f_high); % 返回整个窗口滤波结果 y y(end); % 只输出当前最新滤波点注意如果每个步长都调用一次FFT计算量是O(N log N)1024点完全能接受。但如果模型用变步长、步长很小、实时性要求高这么写就要小心。MATLAB Function模块默认按解释型执行每次调用有开销对于1 kHz采样率没问题到了100 kHz可能就扛不住。第二种用DSP System Toolbox里的Spectrum Analyzer或Spectrum Filter模块不过区别不大。在线FFT滤波的坑主要在于缓冲区长度和滤波输出的延迟。窗口长度决定了频率分辨率分辨率越细需要的N越大延迟也越大。你需要在“能分辨出噪声频率”和“输出延迟不能影响闭环稳定性”之间找平衡。一般测下来N取1024就够用了延迟在50 Hz信号里大概是512个点的相移如果只是给示波器显示没问题但要是给控制器反馈那必须考虑这个延迟甚至要额外做相位补偿。5.3 从Simulink外部模式直接获取波形如果你正在用Simulink的外部模式External mode跑硬件在环或者实时仿真那么数据可以直接从目标机上的信号流里导出。这时FFT滤波可以放在主机端作为后处理也可以作为控制回路里的一个子模块。建议优先用离线方式验证滤波器参数别直接在实时环境里调参数。我见过有人在外部模式下把FFT滤波模块接进控制器反馈通道结果因为窗口延迟导致系统振荡折腾了一下午。原因是1024点的FFT意味着大约1秒的延迟对于一个快速回路是完全不能接受的。6. 常见问题与排查技巧实录这部分才是真正的干货所在。我在做这类项目时几乎每个问题都亲眼见过网上很多帖子也反复出现。整理成速查表方便你踩坑时翻。6.1 频谱泄漏导致“多出来的峰”症状明明信号只有一个50 Hz正弦幅值谱却像一颗彗星周围一堆小峰截止滤波后波形边缘呈波浪状。原因数据点数N不是信号周期的整数倍导致做FFT时隐式周期性延拓把端点强行接续产生不连续跳变。矩形窗的旁瓣大能量向相邻频率漏出。对策优先加汉宁窗或布莱克曼窗减少泄漏。如果必须保留幅值精度用平顶窗flattop。调整N尽量让长度包含整数个信号周期但前提是你知道信号基频。考虑用Zoom-FFT或者Goertzel算法精确计算某个频率而不是完全依赖长FFT。6.2 滤波后波形在突变位置产生“振铃”症状滤波后的矩形波或阶跃波形在跳变沿前后出现明显的高频抖动振荡幅度还不小。原因频域硬截止相当于给频谱乘了一个方窗时域上等效于和sinc函数卷积sinc函数的旁瓣引起过冲振荡也就是吉布斯现象。对策不要用硬截止改用余弦过渡带。比如在截止频率附近构造一个宽度为过渡带宽的余弦斜坡掩码让频谱平滑衰减到零。过渡带宽一般取信号最高有效频率的1/10到1/20比如保留50 Hz过渡带可设5 Hz左右。接受一定的振铃毕竟这是线性时不变滤波的固有代价。如果必须完全避免那就考虑用FIR滤波器设计比如fir1通过窗函数设计实现更平滑的频响。6.3 做完IFFT得到的信号带虚部或共轭不对称症状ifft(Y_filtered)输出是复数real()取实部后发现波形幅值比原始小或者形状失真。原因频谱掩码操作时把正频率置零了但负频率没对称置零破坏共轭对称。也可能是频率轴设计错误把0频位置当成中心处理。对策始终使用完整频率轴f (0:N-1)*(Fs/N)做掩码不要自己单独构造单边频率轴来操作。如果没有把握干脆用fftshift把零频移到中心再操作最后ifftshift回去再IFFT。这个套路不容易错。Y fftshift(fft(x)); f ( -ceil((N-1)/2):floor((N-1)/2) ) * (Fs/N); % 这里的f长度和Y一致 % 用逻辑索引对Y操作 Y_filt Y .* mask; y real(ifft(ifftshift(Y_filt)));6.4 导入的.mat文件变量名不确定代码总是报错症状明明load成功了但脚本里写死了data.x结果文件里变量叫data.signal1运行直接报错找不到字段。对策先用whos(-file, filename)列出变量名再用动态字段访问或递归提取。如果文件里只有一个变量且是结构体可以用S load(filename); fn fieldnames(S); if numel(fn) 1 data_struct S.(fn{1}); end如果结构体内部还有嵌套字段不确定信号在哪个层次写一个小递归遍历函数把所有包含数值数组的字段都找出来然后根据长度和时间轴筛选。我通常会写一个通用findSignal函数参数是结构体和一个最小长度阈值遍历所有数字字段返回最大的那个向量。6.5 采样时间不等距导致频谱异常症状FFT画出来的频谱像经济危机后的股票曲线低频能量特别大明显不合理。幅值谱里0 Hz附近巨大峰但时域明明没有直流。原因Simulink变步长仿真或者硬件采集抖动导致时间轴不均匀。FFT默认等间隔如果原始时间不是等间隔FFT结果完全失去物理意义。对策先检查时间步长dt diff(t)是否大致恒定。方法if std(dt) / mean(dt) 0.01 % 重采样 t_uniform linspace(t(1), t(end), N_target); x_uniform interp1(t, x, t_uniform, spline); end用spline插值在时间点较少但曲线光滑时表现很好点数太多的话用linear更快。重采样后频率分辨率会变但能保证频谱可解释。6.6 FFT计算时间过长卡顿明显症状处理10万点数据直接fft还好但每次循环里都调用就明显慢。对策一是在可能的情况下用2的幂次点数fft最快二是预分配输出数组避免动态扩展三是改到coder.extrinsic调用MATLAB函数前注意离线优化四是如果只是看频谱可以用pwelch它做的是分段平均周期图不是直接对整段FFT性能更好而且更稳。另外提个提醒很多人纠结FFT点数N越大越好。实际上频率分辨率df Fs/N点数翻倍分辨率才提高一倍但计算量增长接近线性对数。对于大多数工程信号N取2048或4096足够看到细节。如果你非要0.01 Hz的分辨率那意味着N至少是Fs/0.01容易变得很慢这时可以考虑先降采样或者分段处理。6.7 滤波器参数如何快速确定我以前总爱问“普通信号用什么截止频率”。实话是没有万能答案但有一个可复现的调试套路画出原始信号的幅值谱。用鼠标取点或者ginput在谱图上点选要保留的频率区间。对选择的区间做掩码滤波后画时域图对比。如果发现某噪声没滤干净看它的频率在哪再缩小保留区间或增加过渡带。如果发现有用信号被削平看保留区间是否太窄扩大。这套“看谱调滤波”的方法比设计滤波器再扫频要快得多尤其在噪声是离散窄带分量时。我就是靠这招从项目中省下大量仿真时间。7. 一个完整示例将Simulink示波器数据和外部.mat数据统一滤波为了让你真正能照着做我把整个流程串成一个实际案例。假设你有一个Simulink模型servo_ctrl.slx模型里有你想分析的输出信号theta另外你手头还有一个外部采集的.mat文件sensor.mat。两者都需要做FFT滤波滤除30 Hz以上的高频抖动保留低频控制信号。第一步从Simulink导出数据model servo_ctrl; load_system(model); out sim(model, StopTime, 10); % 假定模型里存在To Workspace模块输出变量名是 theta_out theta_sim out.theta_out.Data; % 列向量 t_sim out.theta_out.Time; Fs_sim 1 / mean(diff(t_sim));第二步加载外部.mat数据S load(sensor.mat); % 假设里面字段是 measurement结构体里有 t 和 y t_ext S.measurement.t; y_ext S.measurement.y; Fs_ext 1 / median(diff(t_ext));第三步写一个通用滤波函数两者复用function x_filt lowpass_fft(x, Fs, cutoff, transition) N length(x); x x(:) - mean(x); Y fftshift(fft(x)); f ( -ceil((N-1)/2):floor((N-1)/2) ) * (Fs/N); % 构造平滑掩码保留 [-cutoff cutoff]过渡带 transition mask ones(N,1); % 用 sigmoid 或余弦过渡均可 edge cutoff transition/2; inside abs(f) (cutoff - transition/2); outside abs(f) edge; mask(inside) 1; mask(outside) 0; % 过渡带内用线性插值也可以用三次 slopeIdx ~inside ~outside; f_slope f(slopeIdx); % 线性衰减 mask(slopeIdx) (edge - abs(f_slope)) / transition; % 大于edge的部分为0正常 Y_filt Y .* mask; x_filt real(ifft(ifftshift(Y_filt))); % 加窗补偿没有加窗就不需要但去直流会降低整体均值 end第四步分别滤波并绘制对比theta_filt lowpass_fft(theta_sim, Fs_sim, 30, 5); y_filt lowpass_fft(y_ext, Fs_ext, 30, 5); figure; subplot(2,2,1); plot(t_sim, theta_sim); title(Simulink原始); subplot(2,2,2); plot(t_sim, theta_filt); title(Simulink滤波); subplot(2,2,3); plot(t_ext, y_ext); title(外部原始); subplot(2,2,4); plot(t_ext, y_filt); title(外部滤波);这个例子我实测跑过关键点有两个。一是Simulink导出的时间轴t_sim可能不是从0开始的但你做FFT时完全不关心绝对时间起点只要等间隔就行。二是外部.mat里的时间轴如果有毛刺比如丢失了几帧必须先重采样lowpass_fft里没有重采样所以你在调用前一定要保证数据等间隔。这点我建议写成一个检查函数确保万无一失。8. 进阶用零相位滤波替代FFT硬截止时的思考有些场景下你会发现FFT滤波虽然灵活但处理完的波形总归有一点相位延迟如果你用滑动窗口在线滤波延迟更明显。这时候可以考虑MATLAB的filtfilt函数做零相位FIR滤波它是把信号正向和反向各滤波一遍抵消相位偏移。但要注意filtfilt并不是频率域掩码它需要你设计一个滤波器比如巴特沃斯。我在实际项目中的取舍是如果我只是想看看哪些频率成分存在或者要做批量离线分析优先用FFT滤波。如果滤波后的信号要进入控制系统闭环或者后续要对波形做精确对比比如测量相位差我会用filtfilt配合FIR滤波器。如果两者都想兼顾我的最终方案是在FFT域先看频谱确认保留频率区间然后用fir1设计一个同宽带的线性相位FIR滤波器再用filtfilt处理。这样既得到了频谱可视化的指导又获得了零相位、平滑的时域输出。这里提醒一下fir1设计时阶数要选好一般n 3 * (Fs / cutoff)左右太低会过度带太宽太高计算量大。设计完可以用freqz看频响再决定是否调整。9. 性能优化与批处理小技巧如果你的数据特别多比如几百个.mat文件每个都要滤波并保存结果那写一个批处理循环是必然的。这里分享几个能明显提升效率的经验。第一使用dir和正则表达式筛选文件。dir(*.mat)返回结构体数组循环处理时用{files.name}获取文件名再用regexp筛选你要的编号。第二尽量复用FFT的索引和掩码。如果所有数据长度相同、采样率相同那么频率轴和掩码向量可以只算一次放进循环外不要在每次循环里重新生成。第三利用MATLAB的向量化操作。在循环里不要逐点操作尤其是频域掩码直接矩阵乘法或逻辑索引性能差距很大。一次处理一万个点几乎感觉不到卡顿。第四考虑用parfor并行循环。如果你的机器有多核而且各文件之间没有依赖关系直接用parfor替代for能快好几倍。前提是每个迭代中写入独立文件互不干扰。10. 日常经验的一些碎碎念做了这么多FFT滤波我最大的感受是FFT滤波不是魔法它只是把“滤波”问题变成“看频谱做掩码”的问题。它的价值在于可视化、可交互、可试错而它的代价是块处理带来的延迟和边界效应。很多人喜欢写一个“万能滤波函数”到处套结果换一个数据就失真。我的建议是先花十分钟把原始数据的频谱看明白再去选滤波方式这永远是效率最高的路径。另外MATLAB的FFT实现本身就很成熟你不需要去实现一个FFT函数重点放在数据清洗和频谱解释上。遇到波形奇怪先怀疑数据和时间轴再怀疑掩码构造最后才考虑是不是FFT本身的问题。我踩过一次坑明明滤波算法一点问题都没有但原始数据里时间轴有零点漂移导致频谱0 Hz附近异常隆起滤掉低频后波形反而偏离原始形状。那次排查费了我足足半小时后来学乖了所有数据进FFT之前一律先min(t)看看起点和时间差。如果你做的是硬件数据示波器导出的CSV或者.mat文件还要额外注意量化误差和触发噪声。对于非常小幅度的信号量化噪声在频谱上表现为平直的背景底噪FFT滤波对白噪声有一定抑制作用但它无法消除宽带噪声中与信号重叠的部分。真要处理强噪声可以考虑小波变换或自适应滤波那是另一个话题。这个FFT滤波方案还能怎么扩展我目前常做的两个方向一是把滤波后的频谱直接用于自动特征提取比如判断电机轴系是否有故障频率二是结合App Designer写一个小工具让同事直接从示波器数据文件拖进来点击按钮就完成滤波和报告生成省去教别人写脚本的功夫。这些后续有机会再单独展开聊。如果这篇内容能帮你实际解决问题那就已经值了。
返回列表