1. 从零开始为什么选择Scilab处理正弦信号如果你正在信号处理、通信或者生物医学工程领域摸索大概率会频繁遇到一个老朋友正弦信号。无论是模拟一个纯净的音频分析一个振动传感器的输出还是处理心电图ECG这类生物电信号正弦波及其组合都是绕不开的基础。很多朋友上手就直奔MATLAB这当然没问题但今天我想聊聊一个同样强大、却完全开源免费的替代品——Scilab以及我们如何用它来搞定正弦信号处理。我最初接触Scilab是因为一个学生项目预算有限MATLAB的授权费用让我们望而却步。在尝试了PythonNumPy/SciPy和Octave之后我们最终选择了Scilab。原因很简单它的语法和操作逻辑与MATLAB高度相似学习曲线平缓内置了丰富的信号处理工具箱开箱即用最重要的是它的图形界面和绘图能力对于教学和快速原型验证非常友好。处理正弦信号本质上就是和幅度、频率、相位这三个参数打交道进行生成、变换、分析和可视化。Scilab在这些基础操作上提供了极其简洁明了的函数让你能快速把理论公式变成可视化的结果这对于理解概念和调试算法至关重要。所以这篇内容适合谁呢如果你是信号处理的新手想找一个轻量、免费的工具入门或者你是学生、研究人员受限于软件成本亦或是你习惯了MATLAB但需要在一个开源环境中复现或分享你的工作。那么跟着我用Scilab走一遍正弦信号处理的核心流程你会发现实现那些教科书上的经典操作并没有想象中那么复杂。我们将从生成一个正弦波开始一步步深入到频谱分析、滤波应用甚至结合最新的研究热点比如MIMO雷达信号处理中的一些基础仿真思路看看这些基础工具如何支撑更前沿的探索。2. 环境搭建与第一个正弦波从安装到绘图工欲善其事必先利其器。我们的第一步是准备好Scilab这个“工坊”。整个过程非常简单远没有配置一些开源科学计算环境那么繁琐。2.1 Scilab的获取与安装Scilab的官方网站提供了Windows、macOS和Linux系统的安装包。对于Windows和macOS用户直接下载对应的安装程序像安装普通软件一样点击“下一步”即可。Linux用户可以通过包管理器安装例如在Ubuntu上可以使用sudo apt-get install scilab命令。安装完成后打开Scilab你会看到一个类似命令行的主窗口控制台和一个包含文件浏览器、变量查看器等面板的界面。对于信号处理我们绝大部分操作都会在控制台中输入指令或者编写脚本文件.sce 后缀来批量执行。安装后我建议先做一个快速验证在控制台输入11并按回车如果返回ans 2.说明环境运行正常。那个小数点“.”是Scilab默认显示浮点数的方式不必在意。2.2 生成你的第一个正弦信号在信号处理中我们通常在离散的时间点上处理信号。这意味着我们需要定义一个时间向量t它代表了一系列等间隔的采样时刻。然后根据正弦波的公式s(t) A * sin(2*π*f*t φ)来计算每个时刻的信号值。其中A是幅度f是频率单位Hzφ是初始相位单位弧度。假设我们想生成一个幅度为1频率为5 Hz初始相位为0的正弦波采样频率为100 Hz即每秒采集100个点持续时间为1秒。在Scilab控制台中我们可以这样操作// 定义基本参数 A 1; // 幅度 f 5; // 频率 (Hz) Fs 100; // 采样频率 (Hz) T 1; // 持续时间 (秒) phi 0; // 初始相位 (弧度) // 生成时间向量。从0开始到T结束不包含T步长为1/Fs。 // 更常用的方法是使用linspace生成指定数量的点。 N T * Fs; // 总采样点数 t linspace(0, T, N1); // 生成N1个点包含0和T时刻 t t(1:$-1); // 去掉最后一个点得到N个点时间范围为[0, T)这是更常见的做法 // 或者更直接地 t (0:N-1)/Fs; // 生成正弦波信号 s A * sin(2*%pi * f * t phi);这里有几个关键点需要注意。%pi是Scilab内置的圆周率π常数。linspace(0, T, N1)生成了从0到T包含T的N1个等间隔点。但我们通常希望时间范围是[0, T)即包含0但不包含T所以我去掉了最后一个点。另一种更清晰、更高效的做法是直接使用t (0:N-1)/Fs;它明确地表示了第n个采样点对应的时间是n/Fs秒。2.3 可视化让信号“看得见”生成了一堆数字但最直观的方式永远是看图。Scilab的plot函数非常强大。我们可以将刚刚生成的正弦波画出来// 绘制时域波形 clf(); // 清除当前图形窗口 plot(t, s); xgrid(1); // 添加网格线参数1表示启用 xtitle(‘5 Hz正弦波时域图‘ ‘时间 (秒)‘ ‘幅度‘);执行这段代码会弹出一个图形窗口显示出一个标准的正弦波形。你可以清晰地看到它在1秒内完成了5个完整的周期这与我们设定的5 Hz频率相符。通过调整f参数重新生成s并绘图你可以立即观察到频率变化对波形的影响。例如将f改为10你会看到波形更“密集”了。注意在脚本文件中建议始终在绘图前使用clf()来清除图形窗口避免新旧图形叠加在一起造成混淆。xgrid(1)添加的网格线能帮助你更准确地读数。这一步虽然基础但至关重要。它建立了从数学公式到可视化结果的直接桥梁。在后续更复杂的操作中养成“生成-绘图-验证”的习惯能帮你快速定位问题理解算法行为。很多初学者卡在算法实现上往往是因为没有直观地看到中间数据的样子。在Scilab里几乎每一个重要步骤后你都可以用几行绘图代码来检查结果这是它作为教学和原型工具的一大优势。3. 深入正弦信号的核心参数分析与操作当我们能生成一个正弦波后下一步就是“玩弄”它通过改变其参数来观察效果并实现一些基本操作。这不仅是熟悉Scilab函数的过程更是加深对信号本身理解的过程。3.1 幅度、频率与相位的独立控制实验让我们设计一个小实验在一个图上同时绘制三个信号原始信号5 Hz 幅度1 相位0、改变幅度后的信号、改变频率后的信号以及改变相位后的信号。为了清晰我们使用subplot来创建子图。// 基础信号参数 A 1; f 5; Fs 100; T 1; phi 0; N T * Fs; t (0:N-1)/Fs; s_ref A * sin(2*%pi * f * t phi); // 参考信号 // 1. 改变幅度 (A 2) s_amp 2 * sin(2*%pi * f * t phi); // 2. 改变频率 (f 10 Hz) s_freq A * sin(2*%pi * 10 * t phi); // 3. 改变相位 (phi %pi/2 即90度) s_phase A * sin(2*%pi * f * t %pi/2); // 绘制对比图 clf(); subplot(2,2,1); plot(t, s_ref); xtitle(‘参考信号: A1 f5Hz phi0‘ ‘时间(s)‘ ‘幅度‘); xgrid(1); subplot(2,2,2); plot(t, s_amp); xtitle(‘幅度加倍: A2‘ ‘时间(s)‘ ‘幅度‘); xgrid(1); subplot(2,2,3); plot(t, s_freq); xtitle(‘频率加倍: f10Hz‘ ‘时间(s)‘ ‘幅度‘); xgrid(1); subplot(2,2,4); plot(t, s_phase); xtitle(‘相位偏移90度: phiπ/2‘ ‘时间(s)‘ ‘幅度‘); xgrid(1);运行这段代码你会得到四张并列的图。对比它们你可以直观地看到幅度影响波形的“高度”或能量频率影响波形的“疏密”或变化快慢相位影响波形在时间轴上的“起始位置”。理解这三者的独立影响是后续进行信号调制、合成和滤波的基础。3.2 信号的基本运算叠加与调制现实中的信号很少是单一频率的正弦波它们往往是多个正弦波的叠加或者被另一个信号所调制。Scilab中处理这些运算就像做普通数学计算一样简单。信号叠加生成两个不同频率的正弦波然后将它们相加得到一个复合信号。f1 3; f2 20; // 两个频率 s1 sin(2*%pi * f1 * t); s2 0.5 * sin(2*%pi * f2 * t); // 第二个信号幅度小一些 s_sum s1 s2; // 信号叠加 clf(); subplot(3,1,1); plot(t, s1); xtitle(‘3 Hz 信号‘); xgrid(1); subplot(3,1,2); plot(t, s2); xtitle(‘20 Hz 信号‘); xgrid(1); subplot(3,1,3); plot(t, s_sum); xtitle(‘叠加后的信号‘); xgrid(1);从第三张图可以看到叠加后的信号波形变得复杂了它包含了低频3Hz的轮廓和高频20Hz的细节纹路。这就是为什么一段音乐或语音信号看起来如此复杂的原因——它是许多不同频率、不同幅度正弦波的集合。幅度调制AM这是一种简单的调制方式用低频的信号消息去控制高频正弦波载波的幅度。我们可以模拟一个最简单的AM调制。// 载波信号 (高频) fc 50; // 载波频率 50Hz carrier sin(2*%pi * fc * t); // 消息信号 (低频) fm 2; // 消息频率 2Hz message 0.8 * sin(2*%pi * fm * t); // 幅度小于1避免过调制 // AM调制: 已调信号 (1 消息) * 载波 // 1 message 将消息信号偏移到正值区域 s_am (1 message) .* carrier; // 注意是点乘 (.*) clf(); subplot(3,1,1); plot(t, message); xtitle(‘消息信号 (2Hz)‘); xgrid(1); subplot(3,1,2); plot(t, carrier); xtitle(‘载波信号 (50Hz)‘); xgrid(1); subplot(3,1,3); plot(t, s_am); xtitle(‘AM已调信号‘); xgrid(1);注意这里的.*是点乘运算符用于对两个向量进行逐元素相乘。观察AM已调信号你会发现载波正弦波的“包络线”即其幅度变化的轮廓形状与消息信号一致。这就是AM调制的核心信息蕴含在幅度的变化中。实操心得在进行向量运算时务必留意Scilab的运算符。*是矩阵乘法而.*是逐元素乘法。对于同样大小的向量a和ba * b会报错除非是行向量乘列向量做内积而a .* b才是我们通常需要的对应点相乘。这是从MATLAB转过来的用户最容易踩的坑之一。通过这些操作你不仅学会了Scilab的语法更重要的是建立了信号如何通过数学运算进行组合与变换的直观感受。接下来我们将进入更核心的领域如何分析一个信号里到底有哪些频率成分。4. 频谱分析用FFT看清信号的“成分”时域图告诉我们信号随时间如何变化但很多时候我们更关心信号由哪些频率成分构成。例如在音频处理中我们想知道一段声音里高音多还是低音多在故障诊断中我们想从振动信号里找出异常的频率分量。这时就需要傅里叶变换而它的高效算法实现就是快速傅里叶变换FFT。Scilab内置的fft函数让这一切变得轻而易举。4.1 FFT基础与单频信号分析让我们对之前生成的单一频率5 Hz正弦波做FFT看看在频域它是什么样子。// 生成单一频率信号 A 1; f0 5; Fs 100; T 1; N T * Fs; // N100 t (0:N-1)/Fs; s A * sin(2*%pi * f0 * t); // 执行FFT S fft(s); // FFT结果S是复数向量。我们通常关心其幅度谱。 // 取绝对值得到幅度并由于对称性通常只取前一半对于实信号 S_mag abs(S); // 构建对应的频率轴 // 频率分辨率为 Fs/N频率轴从0到Fs实际上到Fs/2即可 freq_axis (0:N-1) * Fs / N; // 绘制双边幅度谱 clf(); subplot(2,1,1); plot(freq_axis, S_mag, ‘-o‘); // 用‘o‘标记数据点 xtitle(‘双边幅度谱‘ ‘频率 (Hz)‘ ‘幅度‘); xgrid(1); // 通常我们更关注0到Fs/2的单边谱 N_half ceil(N/2); // 取前半部分索引 S_mag_single S_mag(1:N_half); freq_axis_single freq_axis(1:N_half); subplot(2,1,2); plot(freq_axis_single, S_mag_single, ‘-o‘); xtitle(‘单边幅度谱‘ ‘频率 (Hz)‘ ‘幅度‘); xgrid(1);运行代码后你会看到在双边谱中幅度在5 Hz和95 Hz即100-5 Hz处有两个尖峰。这是因为对于实值信号其频谱是共轭对称的。单边谱则只显示了0到50 HzFs/2的部分这里只有一个清晰的尖峰在5 Hz处。尖峰的幅度大约是50而不是我们信号时域的幅度1。这是因为FFT的结果没有进行归一化。实际的幅度信息需要将FFT结果除以点数N。对于单频信号归一化后的幅度峰值应为 A/2对于双边谱或 A对于单边谱若将负频率能量合并。我们可以修正一下// 绘制归一化的单边幅度谱 S_mag_normalized S_mag / N; // 双边谱归一化 // 对于单边谱除直流分量外其他频率分量能量是双边谱的两倍所以乘以2 S_mag_single_norm S_mag_normalized(1:N_half); S_mag_single_norm(2:$) 2 * S_mag_single_norm(2:$); // 从第二个点开始乘2第一个是0Hz直流 clf(); plot(freq_axis_single, S_mag_single_norm, ‘-o‘); xtitle(‘归一化单边幅度谱‘ ‘频率 (Hz)‘ ‘幅度‘); xgrid(1); xstring(f0, A0.05, ‘峰值 ≈ ‘ string(A)); // 在峰值处添加文本标注现在你应该能看到在5 Hz处的峰值幅度非常接近1这与我们时域设置的幅度A1相符。这个“归一化”步骤在实际分析中非常重要它能让你从频谱图中直接读出信号成分的真实幅度。4.2 多频信号与频谱泄露现象现在让我们分析一个包含多个频率成分的信号并引入一个关键概念频谱泄露。// 生成一个包含3 Hz和20 Hz的信号 f1 3; f2 20; s_multi sin(2*%pi*f1*t) 0.5*sin(2*%pi*f2*t); // 计算FFT并绘制归一化单边谱 S_multi fft(s_multi); N length(s_multi); S_mag_norm abs(S_multi) / N; freq_axis (0:N-1)*Fs/N; N_half ceil(N/2); S_single_norm S_mag_norm(1:N_half); S_single_norm(2:$) 2 * S_single_norm(2:$); f_single freq_axis(1:N_half); clf(); plot(f_single, S_single_norm, ‘-o‘); xtitle(‘多频信号频谱 (3Hz 20Hz)‘ ‘频率 (Hz)‘ ‘幅度‘); xgrid(1);理想情况下频谱图应该只在3 Hz和20 Hz处有两条干净的竖线。但实际绘图你会发现尖峰底部有“拖尾”能量似乎扩散到了旁边的频率点上。这就是“频谱泄露”。产生的主要原因是我们分析的信号片段1秒不是信号周期的整数倍。3 Hz信号在1秒内恰好有3个完整周期所以它的频谱比较干净。但20 Hz信号在1秒内有20个完整周期也是整数所以泄露也不明显。如果我们把频率改为一个非整数倍周期的值泄露会非常严重。// 使用非整数倍周期的频率 f_leak 5.3; // 5.3 Hz 在1秒内有5.3个周期不是整数 s_leak sin(2*%pi * f_leak * t); // ... (计算并绘制频谱的代码同上)你会发现频谱图上本应在5.3 Hz处的单一尖峰变成了一堆散布在多个频率点上的小山包主峰也不在精确的5.3 Hz上。这严重影响了频率和幅度的估计精度。如何减轻频谱泄露答案是使用窗函数。窗函数在时域对信号两端进行平滑衰减减少因信号截断假设信号在观察窗外周期性重复带来的不连续性。Scilab提供了window函数来生成各种窗。// 使用汉宁窗 (Hanning Window) win window(‘hn‘ N); // ‘hn‘ 代表 Hanning s_windowed s_leak .* win; // 将信号与窗函数点乘 // 分别绘制原始信号和加窗信号的频谱进行对比 // ... (计算两者FFT和归一化单边谱) S_raw fft(s_leak); S_win fft(s_windowed); // 计算幅度谱 (注意加窗后信号能量有损失归一化需考虑窗的相干增益此处仅作对比) S_raw_mag abs(S_raw)/N; S_raw_mag(1:N_half) S_raw_mag(1:N_half) * 2; S_win_mag abs(S_win)/N; S_win_mag(1:N_half) S_win_mag(1:N_half) * 2; clf(); subplot(2,1,1); plot(f_single, S_raw_mag(1:N_half)); xtitle(‘原始信号频谱 (泄露严重)‘ ‘频率 (Hz)‘ ‘幅度‘); xgrid(1); subplot(2,1,2); plot(f_single, S_win_mag(1:N_half)); xtitle(‘加汉宁窗后频谱‘ ‘频率 (Hz)‘ ‘幅度‘); xgrid(1);加窗后你会看到频谱泄露产生的“拖尾”被显著抑制了旁瓣主峰旁边的小峰更低。虽然主峰看起来变“胖”了主瓣宽度增加频率分辨率下降但频率估计更准确幅度估计也受泄露影响更小。这是一种典型的折衷用分辨率换取频谱纯度。在实际工程中选择什么样的窗汉宁窗、汉明窗、布莱克曼窗等取决于你对主瓣宽度和旁瓣衰减的具体要求。实操心得进行FFT分析时务必关注信号长度是否为频率成分周期的整数倍。如果不是一定要考虑加窗。fft函数本身很快但理解其输出结果的含义复数、对称性、归一化、频率轴构建和潜在问题泄露、栅栏效应才是关键。Scilab的fft函数默认不进行任何窗处理把控制权完全交给了用户这既是灵活性的体现也要求使用者具备相关知识。掌握了FFT你就拥有了观察信号频率成分的“显微镜”。接下来我们将利用这个工具对信号进行实际的滤波处理。5. 滤波实战从噪声中提取目标信号在实际应用中信号几乎总是与噪声混杂在一起。滤波就是从混合信号中分离出我们感兴趣部分的过程。Scilab的analog和digital滤波器设计工具箱功能强大但对于入门我们可以从经典的有限长单位冲激响应FIR滤波器开始它易于理解和设计。5.1 设计一个简单的FIR低通滤波器假设我们有一个混合了低频5 Hz和高频30 Hz噪声的信号我们想保留低频部分滤除高频噪声。我们可以设计一个低通滤波器让低于某个截止频率比如15 Hz的信号通过而高于它的信号被衰减。FIR滤波器的一个简单设计方法是使用窗函数法。Scilab的wfir函数可以帮我们完成这个工作。// 设计一个FIR低通滤波器 ftype ‘lp‘; // 低通 forder 64; // 滤波器阶数抽头数-1阶数越高过渡带越陡但延迟和计算量越大 fcut 15/(Fs/2); // 归一化截止频率。Fs/2是奈奎斯特频率。15Hz是我们想要的截止频率。 wtype ‘hn‘; // 使用汉宁窗作为设计窗 hm wfir(ftype, forder, fcut, wtype); // hm是滤波器的冲激响应系数 // 绘制滤波器的频率响应 clf(); [hm_freq, fr] frmag(hm, 512); // 计算频率响应512个点 fr_Hz fr * (Fs/2); // 将归一化频率转换为实际频率(Hz) plot(fr_Hz, 20*log10(hm_freq)); // 纵坐标用分贝(dB)表示 xtitle(‘FIR低通滤波器频率响应 (截止频率15Hz)‘ ‘频率 (Hz)‘ ‘幅度 (dB)‘); xgrid(1); ylim([-80, 5]); // 设置y轴范围方便观察阻带衰减从频率响应图可以看到在15 Hz以下增益接近0 dB即信号基本无衰减通过在15 Hz以上增益迅速下降例如在30 Hz处可能已经衰减了-40 dB或更多这意味着30 Hz的信号幅度会被衰减到原来的1/100以下。这就是低通滤波的效果。5.2 应用滤波器并观察效果现在我们生成一个含噪信号并用设计好的滤波器进行处理。// 生成干净的信号5 Hz正弦波 t (0:999)/Fs; // 生成更长的信号10秒 s_clean sin(2*%pi * 5 * t); // 加入高频噪声30 Hz正弦波作为噪声 noise 0.3 * sin(2*%pi * 30 * t); s_noisy s_clean noise; // 应用滤波器。使用convol函数进行卷积运算。 // 注意卷积会使输出信号变长我们通常取中间部分‘same‘模式或有效部分。 s_filtered convol(hm, s_noisy); // 默认是全卷积输出长度 length(hm)length(s_noisy)-1 // 取中间部分使其长度与输入s_noisy相同 L_hm length(hm); start_idx floor(L_hm/2); s_filtered_cropped s_filtered(start_idx1:start_idxlength(s_noisy)); // 绘制对比图 clf(); subplot(3,1,1); plot(t, s_clean); xtitle(‘原始干净信号 (5 Hz)‘); xgrid(1); subplot(3,1,2); plot(t, s_noisy); xtitle(‘加入30Hz噪声后的信号‘); xgrid(1); subplot(3,1,3); plot(t, s_filtered_cropped); xtitle(‘经过低通滤波后的信号‘); xgrid(1);观察第三幅图你会发现30 Hz的高频纹波噪声几乎被完全去除了波形变得平滑非常接近原始的5 Hz正弦波。这就是滤波器的威力。你也可以通过计算滤波前后信号的频谱来定量观察30 Hz分量是如何被抑制的。注意convol函数进行的是线性卷积。对于FIR滤波更专业的做法是使用filter函数它采用直接II型转置结构更高效且能处理实时流数据。filter函数的用法是y filter(b, 1, x)其中b是滤波器系数向量即我们这里的hm1代表反馈系数向量FIR滤波器没有反馈所以为1。使用filter可以避免手动裁剪信号s_filtered filter(hm, 1, s_noisy); // filter函数输出的信号长度与输入x相同但起始部分存在瞬态响应可以忽略前L_hm个点 s_filtered_valid s_filtered(L_hm1:$);5.3 滤波器的延迟效应细心的你可能发现滤波后的信号波形与原始干净信号在时间上似乎没有完全对齐滤波后的信号好像有轻微的“滞后”。这不是错觉这是FIR滤波器固有的群延迟。对于一个N阶的线性相位FIR滤波器其群延迟是固定的等于(N-1)/2个采样周期。在我们的例子中阶数forder64所以延迟为(64-1)/2 31.5个采样点。在时域图上就表现为波形向右平移。// 计算并补偿群延迟 delay_samples (length(hm)-1)/2; // 滤波器系数的长度是forder1 // 绘制时将滤波后信号的时间轴减去延迟 t_shifted t - delay_samples/Fs; clf(); plot(t, s_clean, ‘b-‘ ‘LineWidth‘ 1.5); // 原始信号蓝色 plot(t_shifted, s_filtered_cropped, ‘r--‘ ‘LineWidth‘ 1.5); // 移位后的滤波信号红色虚线 xtitle(‘滤波信号延迟补偿对比‘ ‘时间 (秒)‘ ‘幅度‘); xgrid(1); legend([‘原始干净信号‘ ‘滤波后信号(时间已补偿)‘]);经过时间补偿后两条曲线应该几乎重合。理解并处理滤波器的延迟在实时处理如音频处理、控制系统中至关重要因为你需要确保处理后的信号与系统其他部分在时间上是同步的。通过这个完整的“设计-应用-分析”流程你不仅学会了如何在Scilab中实现滤波更重要的是理解了滤波器的核心指标截止频率、过渡带、阻带衰减、实际效应噪声抑制、信号延迟以及实现细节卷积 vsfilter函数、群延迟补偿。这些都是将理论应用于实践的关键环节。6. 从仿真到应用连接生物医学与雷达信号处理前沿掌握了正弦信号生成、分析和滤波这些基本功后它们的用武之地远远不止于课堂练习。让我们把视野拓宽看看这些基础工具如何作为基石支撑起像生物医学信号处理和MIMO雷达信号处理这样的前沿领域。6.1 仿真心电图ECG信号与基础处理生物医学信号如心电图ECG本质上是准周期性的生物电信号可以看作是由多个不同频率、不同形态的波形复合而成。我们可以用多个正弦波和特定波形如三角波、高斯波来粗略模拟一个ECG周期并演示基础处理。// 一个简化的ECG周期波形模拟一个心跳 Fs_ecg 360; // ECG常用采样率360 Hz t_heartbeat (0:359)/Fs_ecg; // 模拟1秒的心跳假设心率60bpm // 用几个正弦波分量粗略合成P波、QRS复合波和T波 // 这只是一个非常简化的演示模型 ecg_one_beat 0.5*sin(2*%pi*5*t_heartbeat 0.5) ... // 模拟低频成分 2.0*sin(2*%pi*15*t_heartbeat).*exp(-20*(t_heartbeat-0.2).^2) ... // 模拟QRS波 0.8*sin(2*%pi*3*t_heartbeat - 0.8).*exp(-10*(t_heartbeat-0.6).^2); // 模拟T波 // 生成一段包含多个心跳的ECG信号并加入基线漂移和工频干扰 num_beats 10; ecg_signal []; for i1:num_beats ecg_signal [ecg_signal, ecg_one_beat]; end t_ecg (0:length(ecg_signal)-1)/Fs_ecg; // 加入噪声基线漂移低频和50Hz工频干扰高频 baseline_wander 0.1 * sin(2*%pi*0.2*t_ecg); // 0.2 Hz漂移 powerline_noise 0.05 * sin(2*%pi*50*t_ecg); // 50 Hz干扰 ecg_noisy ecg_signal baseline_wander powerline_noise; // 设计一个带阻滤波器来滤除50Hz工频干扰 // 设计一个窄带阻滤波器陷波滤波器中心频率50Hz f0 50; // 陷波频率 bandwidth 2; // 带宽 (Hz) [b, a] iirnotch(2*f0/Fs_ecg, bandwidth/(Fs_ecg/2)); // iirnotch需要归一化频率 // 应用滤波器 ecg_filtered filter(b, a, ecg_noisy); // 绘图对比 clf(); subplot(3,1,1); plot(t_ecg, ecg_signal); xtitle(‘模拟的干净ECG信号‘); xgrid(1); subplot(3,1,2); plot(t_ecg, ecg_noisy); xtitle(‘加入基线漂移和50Hz干扰的ECG‘); xgrid(1); subplot(3,1,3); plot(t_ecg, ecg_filtered); xtitle(‘经50Hz陷波滤波后的ECG‘); xgrid(1);在这个例子中我们使用了iirnotch函数来设计一个无限长单位冲激响应IIR陷波滤波器专门用于滤除特定频率如50Hz或60Hz工频干扰。与FIR滤波器相比IIR滤波器可以用较低的阶数实现非常尖锐的频率响应但需要注意其相位非线性问题。在生物医学信号处理中滤除工频干扰是预处理中非常常见且关键的一步。通过这个简单的仿真你就能理解实际ECG设备中数字滤波模块所做的事情。6.2 MIMO雷达信号处理中的基础仿真概念多输入多输出MIMO雷达是当前雷达领域的前沿它通过多个发射和接收天线来提升分辨率、抗干扰能力和目标识别能力。其信号处理核心之一就是处理多个发射信号通常是正交的调制信号在目标处反射后在多个接收通道上产生的混合信号。虽然完整的MIMO处理链非常复杂但其基础建模依然离不开正弦信号或更一般的复指数信号和阵列处理。我们可以构建一个极度简化的场景一个具有两个发射天线、两个接收天线的MIMO雷达发射两个频率略有差别的正弦连续波CW来模拟通过频率区分发射通道。// 简化MIMO雷达仿真参数 c 3e8; // 光速 fc 24e9; // 载波频率 24 GHz (典型车载雷达频段) delta_f 10e3; // 发射信号频率差 10 kHz Fs_mimo 1e6; // 采样率 1 MHz T_chirp 1e-3; // 发射脉冲时间 1 ms t_chirp (0:Fs_mimo*T_chirp-1)/Fs_mimo; // 两个发射信号频率分别为 fc 和 fcdelta_f Tx1_signal cos(2*%pi * fc * t_chirp); Tx2_signal cos(2*%pi * (fc delta_f) * t_chirp); // 假设一个目标其回波延迟为 tau R 100; // 目标距离 100米 tau 2*R/c; // 双程延迟 // 两个接收天线接收到的信号忽略幅度衰减、天线方向图等 // 简化模型每个接收信号是两个发射信号延迟后的叠加 Rx1_signal cos(2*%pi * fc * (t_chirp - tau)) cos(2*%pi * (fc delta_f) * (t_chirp - tau)); Rx2_signal cos(2*%pi * fc * (t_chirp - tau)) cos(2*%pi * (fc delta_f) * (t_chirp - tau)); // 假设与Rx1相同实际会有相位差 // 在接收端通过数字下变频和低通滤波得到基带信号 // 1. 数字下变频与发射载频混频 I1 Rx1_signal .* cos(2*%pi * fc * t_chirp); Q1 Rx1_signal .* (-sin(2*%pi * fc * t_chirp)); // 正交分量 // 2. 低通滤波简化此处仅示意性使用移动平均 LPF_len 50; b_lpf ones(1, LPF_len)/LPF_len; // 简单的移动平均滤波器 I1_filtered filter(b_lpf, 1, I1); Q1_filtered filter(b_lpf, 1, Q1); // 基带复信号 BB_signal I1_filtered %i * Q1_filtered; // 对基带信号做FFT观察频率成分 N_fft 2^nextpow2(length(BB_signal)); BB_spectrum fft(BB_signal, N_fft); f_axis (0:N_fft-1)*Fs_mimo/N_fft; // 寻找峰值频率该频率与目标距离有关在FMCW中更直接此处CW仅示意 [mag, idx] max(abs(BB_spectrum(1:N_fft/2))); f_peak f_axis(idx); disp(‘检测到的主要基带频率Hz:‘ f_peak);这个仿真极度简化忽略了MIMO雷达中真正的正交编码、空间角度估计、 Doppler处理等核心环节。但它揭示了一个基本思想复杂的系统级仿真其底层构建模块仍然是正弦信号的生成、调制、混频、滤波和频谱分析。在MIMO雷达中每个发射通道的信号设计、接收通道的回波建模、下变频、以及后续的联合频谱分析如2D-FFT用于距离-速度估计和阵列信号处理如波束成形、DOA估计都可以在Scilab这样的环境中利用我们前面练习过的这些基础操作一步步搭建和验证算法原型。个人体会无论是处理微弱的生物电位还是解析雷达回波信号处理的底层逻辑是相通的。在Scilab中从简单的正弦波开始练习熟练生成、变换、分析和滤波操作就像是练好了扎马步和基本拳法。当面对像ECG去噪或MIMO雷达仿真这类复杂课题时你不会被庞大的系统吓倒因为你知道它们都是由这些基础步骤组合、迭代而成的。你可以快速构建一个简化模型来验证想法这比一开始就陷入复杂的理论公式或庞大的商业仿真软件中要高效得多。我的建议是在掌握了本篇介绍的所有基础操作后尝试用Scilab去复现你专业领域内一篇论文中的某个核心算法框图哪怕只是一个简化版本这个过程会让你对理论和工具的理解产生质的飞跃。 Scilab信号处理入门：从正弦波生成到FFT频谱分析与滤波实战
1. 从零开始为什么选择Scilab处理正弦信号如果你正在信号处理、通信或者生物医学工程领域摸索大概率会频繁遇到一个老朋友正弦信号。无论是模拟一个纯净的音频分析一个振动传感器的输出还是处理心电图ECG这类生物电信号正弦波及其组合都是绕不开的基础。很多朋友上手就直奔MATLAB这当然没问题但今天我想聊聊一个同样强大、却完全开源免费的替代品——Scilab以及我们如何用它来搞定正弦信号处理。我最初接触Scilab是因为一个学生项目预算有限MATLAB的授权费用让我们望而却步。在尝试了PythonNumPy/SciPy和Octave之后我们最终选择了Scilab。原因很简单它的语法和操作逻辑与MATLAB高度相似学习曲线平缓内置了丰富的信号处理工具箱开箱即用最重要的是它的图形界面和绘图能力对于教学和快速原型验证非常友好。处理正弦信号本质上就是和幅度、频率、相位这三个参数打交道进行生成、变换、分析和可视化。Scilab在这些基础操作上提供了极其简洁明了的函数让你能快速把理论公式变成可视化的结果这对于理解概念和调试算法至关重要。所以这篇内容适合谁呢如果你是信号处理的新手想找一个轻量、免费的工具入门或者你是学生、研究人员受限于软件成本亦或是你习惯了MATLAB但需要在一个开源环境中复现或分享你的工作。那么跟着我用Scilab走一遍正弦信号处理的核心流程你会发现实现那些教科书上的经典操作并没有想象中那么复杂。我们将从生成一个正弦波开始一步步深入到频谱分析、滤波应用甚至结合最新的研究热点比如MIMO雷达信号处理中的一些基础仿真思路看看这些基础工具如何支撑更前沿的探索。2. 环境搭建与第一个正弦波从安装到绘图工欲善其事必先利其器。我们的第一步是准备好Scilab这个“工坊”。整个过程非常简单远没有配置一些开源科学计算环境那么繁琐。2.1 Scilab的获取与安装Scilab的官方网站提供了Windows、macOS和Linux系统的安装包。对于Windows和macOS用户直接下载对应的安装程序像安装普通软件一样点击“下一步”即可。Linux用户可以通过包管理器安装例如在Ubuntu上可以使用sudo apt-get install scilab命令。安装完成后打开Scilab你会看到一个类似命令行的主窗口控制台和一个包含文件浏览器、变量查看器等面板的界面。对于信号处理我们绝大部分操作都会在控制台中输入指令或者编写脚本文件.sce 后缀来批量执行。安装后我建议先做一个快速验证在控制台输入11并按回车如果返回ans 2.说明环境运行正常。那个小数点“.”是Scilab默认显示浮点数的方式不必在意。2.2 生成你的第一个正弦信号在信号处理中我们通常在离散的时间点上处理信号。这意味着我们需要定义一个时间向量t它代表了一系列等间隔的采样时刻。然后根据正弦波的公式s(t) A * sin(2*π*f*t φ)来计算每个时刻的信号值。其中A是幅度f是频率单位Hzφ是初始相位单位弧度。假设我们想生成一个幅度为1频率为5 Hz初始相位为0的正弦波采样频率为100 Hz即每秒采集100个点持续时间为1秒。在Scilab控制台中我们可以这样操作// 定义基本参数 A 1; // 幅度 f 5; // 频率 (Hz) Fs 100; // 采样频率 (Hz) T 1; // 持续时间 (秒) phi 0; // 初始相位 (弧度) // 生成时间向量。从0开始到T结束不包含T步长为1/Fs。 // 更常用的方法是使用linspace生成指定数量的点。 N T * Fs; // 总采样点数 t linspace(0, T, N1); // 生成N1个点包含0和T时刻 t t(1:$-1); // 去掉最后一个点得到N个点时间范围为[0, T)这是更常见的做法 // 或者更直接地 t (0:N-1)/Fs; // 生成正弦波信号 s A * sin(2*%pi * f * t phi);这里有几个关键点需要注意。%pi是Scilab内置的圆周率π常数。linspace(0, T, N1)生成了从0到T包含T的N1个等间隔点。但我们通常希望时间范围是[0, T)即包含0但不包含T所以我去掉了最后一个点。另一种更清晰、更高效的做法是直接使用t (0:N-1)/Fs;它明确地表示了第n个采样点对应的时间是n/Fs秒。2.3 可视化让信号“看得见”生成了一堆数字但最直观的方式永远是看图。Scilab的plot函数非常强大。我们可以将刚刚生成的正弦波画出来// 绘制时域波形 clf(); // 清除当前图形窗口 plot(t, s); xgrid(1); // 添加网格线参数1表示启用 xtitle(‘5 Hz正弦波时域图‘ ‘时间 (秒)‘ ‘幅度‘);执行这段代码会弹出一个图形窗口显示出一个标准的正弦波形。你可以清晰地看到它在1秒内完成了5个完整的周期这与我们设定的5 Hz频率相符。通过调整f参数重新生成s并绘图你可以立即观察到频率变化对波形的影响。例如将f改为10你会看到波形更“密集”了。注意在脚本文件中建议始终在绘图前使用clf()来清除图形窗口避免新旧图形叠加在一起造成混淆。xgrid(1)添加的网格线能帮助你更准确地读数。这一步虽然基础但至关重要。它建立了从数学公式到可视化结果的直接桥梁。在后续更复杂的操作中养成“生成-绘图-验证”的习惯能帮你快速定位问题理解算法行为。很多初学者卡在算法实现上往往是因为没有直观地看到中间数据的样子。在Scilab里几乎每一个重要步骤后你都可以用几行绘图代码来检查结果这是它作为教学和原型工具的一大优势。3. 深入正弦信号的核心参数分析与操作当我们能生成一个正弦波后下一步就是“玩弄”它通过改变其参数来观察效果并实现一些基本操作。这不仅是熟悉Scilab函数的过程更是加深对信号本身理解的过程。3.1 幅度、频率与相位的独立控制实验让我们设计一个小实验在一个图上同时绘制三个信号原始信号5 Hz 幅度1 相位0、改变幅度后的信号、改变频率后的信号以及改变相位后的信号。为了清晰我们使用subplot来创建子图。// 基础信号参数 A 1; f 5; Fs 100; T 1; phi 0; N T * Fs; t (0:N-1)/Fs; s_ref A * sin(2*%pi * f * t phi); // 参考信号 // 1. 改变幅度 (A 2) s_amp 2 * sin(2*%pi * f * t phi); // 2. 改变频率 (f 10 Hz) s_freq A * sin(2*%pi * 10 * t phi); // 3. 改变相位 (phi %pi/2 即90度) s_phase A * sin(2*%pi * f * t %pi/2); // 绘制对比图 clf(); subplot(2,2,1); plot(t, s_ref); xtitle(‘参考信号: A1 f5Hz phi0‘ ‘时间(s)‘ ‘幅度‘); xgrid(1); subplot(2,2,2); plot(t, s_amp); xtitle(‘幅度加倍: A2‘ ‘时间(s)‘ ‘幅度‘); xgrid(1); subplot(2,2,3); plot(t, s_freq); xtitle(‘频率加倍: f10Hz‘ ‘时间(s)‘ ‘幅度‘); xgrid(1); subplot(2,2,4); plot(t, s_phase); xtitle(‘相位偏移90度: phiπ/2‘ ‘时间(s)‘ ‘幅度‘); xgrid(1);运行这段代码你会得到四张并列的图。对比它们你可以直观地看到幅度影响波形的“高度”或能量频率影响波形的“疏密”或变化快慢相位影响波形在时间轴上的“起始位置”。理解这三者的独立影响是后续进行信号调制、合成和滤波的基础。3.2 信号的基本运算叠加与调制现实中的信号很少是单一频率的正弦波它们往往是多个正弦波的叠加或者被另一个信号所调制。Scilab中处理这些运算就像做普通数学计算一样简单。信号叠加生成两个不同频率的正弦波然后将它们相加得到一个复合信号。f1 3; f2 20; // 两个频率 s1 sin(2*%pi * f1 * t); s2 0.5 * sin(2*%pi * f2 * t); // 第二个信号幅度小一些 s_sum s1 s2; // 信号叠加 clf(); subplot(3,1,1); plot(t, s1); xtitle(‘3 Hz 信号‘); xgrid(1); subplot(3,1,2); plot(t, s2); xtitle(‘20 Hz 信号‘); xgrid(1); subplot(3,1,3); plot(t, s_sum); xtitle(‘叠加后的信号‘); xgrid(1);从第三张图可以看到叠加后的信号波形变得复杂了它包含了低频3Hz的轮廓和高频20Hz的细节纹路。这就是为什么一段音乐或语音信号看起来如此复杂的原因——它是许多不同频率、不同幅度正弦波的集合。幅度调制AM这是一种简单的调制方式用低频的信号消息去控制高频正弦波载波的幅度。我们可以模拟一个最简单的AM调制。// 载波信号 (高频) fc 50; // 载波频率 50Hz carrier sin(2*%pi * fc * t); // 消息信号 (低频) fm 2; // 消息频率 2Hz message 0.8 * sin(2*%pi * fm * t); // 幅度小于1避免过调制 // AM调制: 已调信号 (1 消息) * 载波 // 1 message 将消息信号偏移到正值区域 s_am (1 message) .* carrier; // 注意是点乘 (.*) clf(); subplot(3,1,1); plot(t, message); xtitle(‘消息信号 (2Hz)‘); xgrid(1); subplot(3,1,2); plot(t, carrier); xtitle(‘载波信号 (50Hz)‘); xgrid(1); subplot(3,1,3); plot(t, s_am); xtitle(‘AM已调信号‘); xgrid(1);注意这里的.*是点乘运算符用于对两个向量进行逐元素相乘。观察AM已调信号你会发现载波正弦波的“包络线”即其幅度变化的轮廓形状与消息信号一致。这就是AM调制的核心信息蕴含在幅度的变化中。实操心得在进行向量运算时务必留意Scilab的运算符。*是矩阵乘法而.*是逐元素乘法。对于同样大小的向量a和ba * b会报错除非是行向量乘列向量做内积而a .* b才是我们通常需要的对应点相乘。这是从MATLAB转过来的用户最容易踩的坑之一。通过这些操作你不仅学会了Scilab的语法更重要的是建立了信号如何通过数学运算进行组合与变换的直观感受。接下来我们将进入更核心的领域如何分析一个信号里到底有哪些频率成分。4. 频谱分析用FFT看清信号的“成分”时域图告诉我们信号随时间如何变化但很多时候我们更关心信号由哪些频率成分构成。例如在音频处理中我们想知道一段声音里高音多还是低音多在故障诊断中我们想从振动信号里找出异常的频率分量。这时就需要傅里叶变换而它的高效算法实现就是快速傅里叶变换FFT。Scilab内置的fft函数让这一切变得轻而易举。4.1 FFT基础与单频信号分析让我们对之前生成的单一频率5 Hz正弦波做FFT看看在频域它是什么样子。// 生成单一频率信号 A 1; f0 5; Fs 100; T 1; N T * Fs; // N100 t (0:N-1)/Fs; s A * sin(2*%pi * f0 * t); // 执行FFT S fft(s); // FFT结果S是复数向量。我们通常关心其幅度谱。 // 取绝对值得到幅度并由于对称性通常只取前一半对于实信号 S_mag abs(S); // 构建对应的频率轴 // 频率分辨率为 Fs/N频率轴从0到Fs实际上到Fs/2即可 freq_axis (0:N-1) * Fs / N; // 绘制双边幅度谱 clf(); subplot(2,1,1); plot(freq_axis, S_mag, ‘-o‘); // 用‘o‘标记数据点 xtitle(‘双边幅度谱‘ ‘频率 (Hz)‘ ‘幅度‘); xgrid(1); // 通常我们更关注0到Fs/2的单边谱 N_half ceil(N/2); // 取前半部分索引 S_mag_single S_mag(1:N_half); freq_axis_single freq_axis(1:N_half); subplot(2,1,2); plot(freq_axis_single, S_mag_single, ‘-o‘); xtitle(‘单边幅度谱‘ ‘频率 (Hz)‘ ‘幅度‘); xgrid(1);运行代码后你会看到在双边谱中幅度在5 Hz和95 Hz即100-5 Hz处有两个尖峰。这是因为对于实值信号其频谱是共轭对称的。单边谱则只显示了0到50 HzFs/2的部分这里只有一个清晰的尖峰在5 Hz处。尖峰的幅度大约是50而不是我们信号时域的幅度1。这是因为FFT的结果没有进行归一化。实际的幅度信息需要将FFT结果除以点数N。对于单频信号归一化后的幅度峰值应为 A/2对于双边谱或 A对于单边谱若将负频率能量合并。我们可以修正一下// 绘制归一化的单边幅度谱 S_mag_normalized S_mag / N; // 双边谱归一化 // 对于单边谱除直流分量外其他频率分量能量是双边谱的两倍所以乘以2 S_mag_single_norm S_mag_normalized(1:N_half); S_mag_single_norm(2:$) 2 * S_mag_single_norm(2:$); // 从第二个点开始乘2第一个是0Hz直流 clf(); plot(freq_axis_single, S_mag_single_norm, ‘-o‘); xtitle(‘归一化单边幅度谱‘ ‘频率 (Hz)‘ ‘幅度‘); xgrid(1); xstring(f0, A0.05, ‘峰值 ≈ ‘ string(A)); // 在峰值处添加文本标注现在你应该能看到在5 Hz处的峰值幅度非常接近1这与我们时域设置的幅度A1相符。这个“归一化”步骤在实际分析中非常重要它能让你从频谱图中直接读出信号成分的真实幅度。4.2 多频信号与频谱泄露现象现在让我们分析一个包含多个频率成分的信号并引入一个关键概念频谱泄露。// 生成一个包含3 Hz和20 Hz的信号 f1 3; f2 20; s_multi sin(2*%pi*f1*t) 0.5*sin(2*%pi*f2*t); // 计算FFT并绘制归一化单边谱 S_multi fft(s_multi); N length(s_multi); S_mag_norm abs(S_multi) / N; freq_axis (0:N-1)*Fs/N; N_half ceil(N/2); S_single_norm S_mag_norm(1:N_half); S_single_norm(2:$) 2 * S_single_norm(2:$); f_single freq_axis(1:N_half); clf(); plot(f_single, S_single_norm, ‘-o‘); xtitle(‘多频信号频谱 (3Hz 20Hz)‘ ‘频率 (Hz)‘ ‘幅度‘); xgrid(1);理想情况下频谱图应该只在3 Hz和20 Hz处有两条干净的竖线。但实际绘图你会发现尖峰底部有“拖尾”能量似乎扩散到了旁边的频率点上。这就是“频谱泄露”。产生的主要原因是我们分析的信号片段1秒不是信号周期的整数倍。3 Hz信号在1秒内恰好有3个完整周期所以它的频谱比较干净。但20 Hz信号在1秒内有20个完整周期也是整数所以泄露也不明显。如果我们把频率改为一个非整数倍周期的值泄露会非常严重。// 使用非整数倍周期的频率 f_leak 5.3; // 5.3 Hz 在1秒内有5.3个周期不是整数 s_leak sin(2*%pi * f_leak * t); // ... (计算并绘制频谱的代码同上)你会发现频谱图上本应在5.3 Hz处的单一尖峰变成了一堆散布在多个频率点上的小山包主峰也不在精确的5.3 Hz上。这严重影响了频率和幅度的估计精度。如何减轻频谱泄露答案是使用窗函数。窗函数在时域对信号两端进行平滑衰减减少因信号截断假设信号在观察窗外周期性重复带来的不连续性。Scilab提供了window函数来生成各种窗。// 使用汉宁窗 (Hanning Window) win window(‘hn‘ N); // ‘hn‘ 代表 Hanning s_windowed s_leak .* win; // 将信号与窗函数点乘 // 分别绘制原始信号和加窗信号的频谱进行对比 // ... (计算两者FFT和归一化单边谱) S_raw fft(s_leak); S_win fft(s_windowed); // 计算幅度谱 (注意加窗后信号能量有损失归一化需考虑窗的相干增益此处仅作对比) S_raw_mag abs(S_raw)/N; S_raw_mag(1:N_half) S_raw_mag(1:N_half) * 2; S_win_mag abs(S_win)/N; S_win_mag(1:N_half) S_win_mag(1:N_half) * 2; clf(); subplot(2,1,1); plot(f_single, S_raw_mag(1:N_half)); xtitle(‘原始信号频谱 (泄露严重)‘ ‘频率 (Hz)‘ ‘幅度‘); xgrid(1); subplot(2,1,2); plot(f_single, S_win_mag(1:N_half)); xtitle(‘加汉宁窗后频谱‘ ‘频率 (Hz)‘ ‘幅度‘); xgrid(1);加窗后你会看到频谱泄露产生的“拖尾”被显著抑制了旁瓣主峰旁边的小峰更低。虽然主峰看起来变“胖”了主瓣宽度增加频率分辨率下降但频率估计更准确幅度估计也受泄露影响更小。这是一种典型的折衷用分辨率换取频谱纯度。在实际工程中选择什么样的窗汉宁窗、汉明窗、布莱克曼窗等取决于你对主瓣宽度和旁瓣衰减的具体要求。实操心得进行FFT分析时务必关注信号长度是否为频率成分周期的整数倍。如果不是一定要考虑加窗。fft函数本身很快但理解其输出结果的含义复数、对称性、归一化、频率轴构建和潜在问题泄露、栅栏效应才是关键。Scilab的fft函数默认不进行任何窗处理把控制权完全交给了用户这既是灵活性的体现也要求使用者具备相关知识。掌握了FFT你就拥有了观察信号频率成分的“显微镜”。接下来我们将利用这个工具对信号进行实际的滤波处理。5. 滤波实战从噪声中提取目标信号在实际应用中信号几乎总是与噪声混杂在一起。滤波就是从混合信号中分离出我们感兴趣部分的过程。Scilab的analog和digital滤波器设计工具箱功能强大但对于入门我们可以从经典的有限长单位冲激响应FIR滤波器开始它易于理解和设计。5.1 设计一个简单的FIR低通滤波器假设我们有一个混合了低频5 Hz和高频30 Hz噪声的信号我们想保留低频部分滤除高频噪声。我们可以设计一个低通滤波器让低于某个截止频率比如15 Hz的信号通过而高于它的信号被衰减。FIR滤波器的一个简单设计方法是使用窗函数法。Scilab的wfir函数可以帮我们完成这个工作。// 设计一个FIR低通滤波器 ftype ‘lp‘; // 低通 forder 64; // 滤波器阶数抽头数-1阶数越高过渡带越陡但延迟和计算量越大 fcut 15/(Fs/2); // 归一化截止频率。Fs/2是奈奎斯特频率。15Hz是我们想要的截止频率。 wtype ‘hn‘; // 使用汉宁窗作为设计窗 hm wfir(ftype, forder, fcut, wtype); // hm是滤波器的冲激响应系数 // 绘制滤波器的频率响应 clf(); [hm_freq, fr] frmag(hm, 512); // 计算频率响应512个点 fr_Hz fr * (Fs/2); // 将归一化频率转换为实际频率(Hz) plot(fr_Hz, 20*log10(hm_freq)); // 纵坐标用分贝(dB)表示 xtitle(‘FIR低通滤波器频率响应 (截止频率15Hz)‘ ‘频率 (Hz)‘ ‘幅度 (dB)‘); xgrid(1); ylim([-80, 5]); // 设置y轴范围方便观察阻带衰减从频率响应图可以看到在15 Hz以下增益接近0 dB即信号基本无衰减通过在15 Hz以上增益迅速下降例如在30 Hz处可能已经衰减了-40 dB或更多这意味着30 Hz的信号幅度会被衰减到原来的1/100以下。这就是低通滤波的效果。5.2 应用滤波器并观察效果现在我们生成一个含噪信号并用设计好的滤波器进行处理。// 生成干净的信号5 Hz正弦波 t (0:999)/Fs; // 生成更长的信号10秒 s_clean sin(2*%pi * 5 * t); // 加入高频噪声30 Hz正弦波作为噪声 noise 0.3 * sin(2*%pi * 30 * t); s_noisy s_clean noise; // 应用滤波器。使用convol函数进行卷积运算。 // 注意卷积会使输出信号变长我们通常取中间部分‘same‘模式或有效部分。 s_filtered convol(hm, s_noisy); // 默认是全卷积输出长度 length(hm)length(s_noisy)-1 // 取中间部分使其长度与输入s_noisy相同 L_hm length(hm); start_idx floor(L_hm/2); s_filtered_cropped s_filtered(start_idx1:start_idxlength(s_noisy)); // 绘制对比图 clf(); subplot(3,1,1); plot(t, s_clean); xtitle(‘原始干净信号 (5 Hz)‘); xgrid(1); subplot(3,1,2); plot(t, s_noisy); xtitle(‘加入30Hz噪声后的信号‘); xgrid(1); subplot(3,1,3); plot(t, s_filtered_cropped); xtitle(‘经过低通滤波后的信号‘); xgrid(1);观察第三幅图你会发现30 Hz的高频纹波噪声几乎被完全去除了波形变得平滑非常接近原始的5 Hz正弦波。这就是滤波器的威力。你也可以通过计算滤波前后信号的频谱来定量观察30 Hz分量是如何被抑制的。注意convol函数进行的是线性卷积。对于FIR滤波更专业的做法是使用filter函数它采用直接II型转置结构更高效且能处理实时流数据。filter函数的用法是y filter(b, 1, x)其中b是滤波器系数向量即我们这里的hm1代表反馈系数向量FIR滤波器没有反馈所以为1。使用filter可以避免手动裁剪信号s_filtered filter(hm, 1, s_noisy); // filter函数输出的信号长度与输入x相同但起始部分存在瞬态响应可以忽略前L_hm个点 s_filtered_valid s_filtered(L_hm1:$);5.3 滤波器的延迟效应细心的你可能发现滤波后的信号波形与原始干净信号在时间上似乎没有完全对齐滤波后的信号好像有轻微的“滞后”。这不是错觉这是FIR滤波器固有的群延迟。对于一个N阶的线性相位FIR滤波器其群延迟是固定的等于(N-1)/2个采样周期。在我们的例子中阶数forder64所以延迟为(64-1)/2 31.5个采样点。在时域图上就表现为波形向右平移。// 计算并补偿群延迟 delay_samples (length(hm)-1)/2; // 滤波器系数的长度是forder1 // 绘制时将滤波后信号的时间轴减去延迟 t_shifted t - delay_samples/Fs; clf(); plot(t, s_clean, ‘b-‘ ‘LineWidth‘ 1.5); // 原始信号蓝色 plot(t_shifted, s_filtered_cropped, ‘r--‘ ‘LineWidth‘ 1.5); // 移位后的滤波信号红色虚线 xtitle(‘滤波信号延迟补偿对比‘ ‘时间 (秒)‘ ‘幅度‘); xgrid(1); legend([‘原始干净信号‘ ‘滤波后信号(时间已补偿)‘]);经过时间补偿后两条曲线应该几乎重合。理解并处理滤波器的延迟在实时处理如音频处理、控制系统中至关重要因为你需要确保处理后的信号与系统其他部分在时间上是同步的。通过这个完整的“设计-应用-分析”流程你不仅学会了如何在Scilab中实现滤波更重要的是理解了滤波器的核心指标截止频率、过渡带、阻带衰减、实际效应噪声抑制、信号延迟以及实现细节卷积 vsfilter函数、群延迟补偿。这些都是将理论应用于实践的关键环节。6. 从仿真到应用连接生物医学与雷达信号处理前沿掌握了正弦信号生成、分析和滤波这些基本功后它们的用武之地远远不止于课堂练习。让我们把视野拓宽看看这些基础工具如何作为基石支撑起像生物医学信号处理和MIMO雷达信号处理这样的前沿领域。6.1 仿真心电图ECG信号与基础处理生物医学信号如心电图ECG本质上是准周期性的生物电信号可以看作是由多个不同频率、不同形态的波形复合而成。我们可以用多个正弦波和特定波形如三角波、高斯波来粗略模拟一个ECG周期并演示基础处理。// 一个简化的ECG周期波形模拟一个心跳 Fs_ecg 360; // ECG常用采样率360 Hz t_heartbeat (0:359)/Fs_ecg; // 模拟1秒的心跳假设心率60bpm // 用几个正弦波分量粗略合成P波、QRS复合波和T波 // 这只是一个非常简化的演示模型 ecg_one_beat 0.5*sin(2*%pi*5*t_heartbeat 0.5) ... // 模拟低频成分 2.0*sin(2*%pi*15*t_heartbeat).*exp(-20*(t_heartbeat-0.2).^2) ... // 模拟QRS波 0.8*sin(2*%pi*3*t_heartbeat - 0.8).*exp(-10*(t_heartbeat-0.6).^2); // 模拟T波 // 生成一段包含多个心跳的ECG信号并加入基线漂移和工频干扰 num_beats 10; ecg_signal []; for i1:num_beats ecg_signal [ecg_signal, ecg_one_beat]; end t_ecg (0:length(ecg_signal)-1)/Fs_ecg; // 加入噪声基线漂移低频和50Hz工频干扰高频 baseline_wander 0.1 * sin(2*%pi*0.2*t_ecg); // 0.2 Hz漂移 powerline_noise 0.05 * sin(2*%pi*50*t_ecg); // 50 Hz干扰 ecg_noisy ecg_signal baseline_wander powerline_noise; // 设计一个带阻滤波器来滤除50Hz工频干扰 // 设计一个窄带阻滤波器陷波滤波器中心频率50Hz f0 50; // 陷波频率 bandwidth 2; // 带宽 (Hz) [b, a] iirnotch(2*f0/Fs_ecg, bandwidth/(Fs_ecg/2)); // iirnotch需要归一化频率 // 应用滤波器 ecg_filtered filter(b, a, ecg_noisy); // 绘图对比 clf(); subplot(3,1,1); plot(t_ecg, ecg_signal); xtitle(‘模拟的干净ECG信号‘); xgrid(1); subplot(3,1,2); plot(t_ecg, ecg_noisy); xtitle(‘加入基线漂移和50Hz干扰的ECG‘); xgrid(1); subplot(3,1,3); plot(t_ecg, ecg_filtered); xtitle(‘经50Hz陷波滤波后的ECG‘); xgrid(1);在这个例子中我们使用了iirnotch函数来设计一个无限长单位冲激响应IIR陷波滤波器专门用于滤除特定频率如50Hz或60Hz工频干扰。与FIR滤波器相比IIR滤波器可以用较低的阶数实现非常尖锐的频率响应但需要注意其相位非线性问题。在生物医学信号处理中滤除工频干扰是预处理中非常常见且关键的一步。通过这个简单的仿真你就能理解实际ECG设备中数字滤波模块所做的事情。6.2 MIMO雷达信号处理中的基础仿真概念多输入多输出MIMO雷达是当前雷达领域的前沿它通过多个发射和接收天线来提升分辨率、抗干扰能力和目标识别能力。其信号处理核心之一就是处理多个发射信号通常是正交的调制信号在目标处反射后在多个接收通道上产生的混合信号。虽然完整的MIMO处理链非常复杂但其基础建模依然离不开正弦信号或更一般的复指数信号和阵列处理。我们可以构建一个极度简化的场景一个具有两个发射天线、两个接收天线的MIMO雷达发射两个频率略有差别的正弦连续波CW来模拟通过频率区分发射通道。// 简化MIMO雷达仿真参数 c 3e8; // 光速 fc 24e9; // 载波频率 24 GHz (典型车载雷达频段) delta_f 10e3; // 发射信号频率差 10 kHz Fs_mimo 1e6; // 采样率 1 MHz T_chirp 1e-3; // 发射脉冲时间 1 ms t_chirp (0:Fs_mimo*T_chirp-1)/Fs_mimo; // 两个发射信号频率分别为 fc 和 fcdelta_f Tx1_signal cos(2*%pi * fc * t_chirp); Tx2_signal cos(2*%pi * (fc delta_f) * t_chirp); // 假设一个目标其回波延迟为 tau R 100; // 目标距离 100米 tau 2*R/c; // 双程延迟 // 两个接收天线接收到的信号忽略幅度衰减、天线方向图等 // 简化模型每个接收信号是两个发射信号延迟后的叠加 Rx1_signal cos(2*%pi * fc * (t_chirp - tau)) cos(2*%pi * (fc delta_f) * (t_chirp - tau)); Rx2_signal cos(2*%pi * fc * (t_chirp - tau)) cos(2*%pi * (fc delta_f) * (t_chirp - tau)); // 假设与Rx1相同实际会有相位差 // 在接收端通过数字下变频和低通滤波得到基带信号 // 1. 数字下变频与发射载频混频 I1 Rx1_signal .* cos(2*%pi * fc * t_chirp); Q1 Rx1_signal .* (-sin(2*%pi * fc * t_chirp)); // 正交分量 // 2. 低通滤波简化此处仅示意性使用移动平均 LPF_len 50; b_lpf ones(1, LPF_len)/LPF_len; // 简单的移动平均滤波器 I1_filtered filter(b_lpf, 1, I1); Q1_filtered filter(b_lpf, 1, Q1); // 基带复信号 BB_signal I1_filtered %i * Q1_filtered; // 对基带信号做FFT观察频率成分 N_fft 2^nextpow2(length(BB_signal)); BB_spectrum fft(BB_signal, N_fft); f_axis (0:N_fft-1)*Fs_mimo/N_fft; // 寻找峰值频率该频率与目标距离有关在FMCW中更直接此处CW仅示意 [mag, idx] max(abs(BB_spectrum(1:N_fft/2))); f_peak f_axis(idx); disp(‘检测到的主要基带频率Hz:‘ f_peak);这个仿真极度简化忽略了MIMO雷达中真正的正交编码、空间角度估计、 Doppler处理等核心环节。但它揭示了一个基本思想复杂的系统级仿真其底层构建模块仍然是正弦信号的生成、调制、混频、滤波和频谱分析。在MIMO雷达中每个发射通道的信号设计、接收通道的回波建模、下变频、以及后续的联合频谱分析如2D-FFT用于距离-速度估计和阵列信号处理如波束成形、DOA估计都可以在Scilab这样的环境中利用我们前面练习过的这些基础操作一步步搭建和验证算法原型。个人体会无论是处理微弱的生物电位还是解析雷达回波信号处理的底层逻辑是相通的。在Scilab中从简单的正弦波开始练习熟练生成、变换、分析和滤波操作就像是练好了扎马步和基本拳法。当面对像ECG去噪或MIMO雷达仿真这类复杂课题时你不会被庞大的系统吓倒因为你知道它们都是由这些基础步骤组合、迭代而成的。你可以快速构建一个简化模型来验证想法这比一开始就陷入复杂的理论公式或庞大的商业仿真软件中要高效得多。我的建议是在掌握了本篇介绍的所有基础操作后尝试用Scilab去复现你专业领域内一篇论文中的某个核心算法框图哪怕只是一个简化版本这个过程会让你对理论和工具的理解产生质的飞跃。