简介围绕实验五“数字信号处理算法实验”这份资源聚焦FIR滤波器、IIR滤波器、FFT及重叠相加法实现FIR滤波四大核心主题面向数字信号处理初学者或正在完成DSP课程实验的学生可支撑滤波器设计、频谱分析和实时滤波等场景的学习与复现。压缩包共217个文件约2.82MB主体为C语言源码与头文件.c/.h及DSP工程构建文件.doj/.dpj/.out另有.dat数据文件、.mak构建脚本、.pcf/.xml配置和.txt说明文本并附一份dsp实验五.docx实验报告结构清晰便于按代码、配置和文档分类查阅。已有790人学习下载。使用这套资料可以系统掌握窗函数法或巴特沃斯滤波器设计、FFT频谱分析以及重叠相加法的实际编程流程源码中包含音频块处理与SRU初始化等模块结合报告中的理论介绍、性能评估与排错思路适合对照实验指导书逐步调试也适合后续移植或扩展DSP项目。1. 一个实验文件里藏着的四个DSP基本功拿到“实验五 数字信号处理算法实验”这套工程文件时第一反应是initSRU.c和blockProcess_audio.c这两类文件重复出现多次像是工程模板里被反复拷贝的骨架。把文件打开看过之后会发现真正值得研究的是四个算法FIR滤波器、IIR滤波器、FFT和重叠相加法。这四个东西恰好覆盖了数字信号处理从频域设计到实时实现的主要路径。FIR解决线性相位和稳定性的问题IIR解决计算效率的问题FFT把频域分析从理论变成工程工具重叠相加法则是把FFT搬回时域做线性卷积的经典手段。这套实验适合两类人一类是刚学完数字信号处理理论、想在C代码里看到滤波器长什么样的人另一类是在嵌入式平台上做音频或传感器处理、需要把算法从MATLAB搬到定点芯片上的工程师。下面先讲FIR因为它是理解后面所有内容的基础。2. FIR滤波器DSP实现从系数设计到C语言定点化2.1 为什么FIR在DSP里这么常用FIR的冲激响应是有限长度的系统输出只依赖当前和过去N-1个输入样本对应差分方程为y[n] b0·x[n] b1·x[n-1] ... b(N-1)·x[n-N1]这里的b0到b(N-1)就是滤波器系数也叫抽头系数。从工程角度看FIR有三个特性让它成为DSP课程的常客。第一它天然稳定因为结构中没有反馈回路极点全在原点第二只要系数具有对称性就能实现线性相位这对音频信号来说意味着群延迟恒定波形不产生相位失真第三滤波过程本质上是做卷积而卷积可以通过FFT加速为后面的重叠相加法埋下伏笔。在blockProcess_audio.c这类实时处理框架里FIR通常会在中断或音频回调中被逐帧调用因此代码层面的循环效率直接决定系统能不能实时跑完。2.2 窗函数法设计滤波器的参数要点设计FIR最常用的是窗函数法。以低通滤波器为例理想低通的单位冲激响应是无限长的sinc函数窗函数法的思路就是截断它。MATLAB里一条命令就能算出系数fs 16000; % 采样率 16kHz fc 3000; % 截止频率 3kHz N 51; % 阶数选择51阶代表51个系数 b fir1(N-1, fc/(fs/2), low, hann(N)); freqz(b, 1, 512, fs); % 查看幅频和相频响应这里N-1是滤波器阶数N是系数个数。fc/(fs/2)是归一化截止频率奈奎斯特频率是8000Hz所以3000Hz对应0.375。窗函数选择hann窗过渡带宽度和阻带衰减是权衡关系窗越长过渡带越窄但计算量变大。如果换成blackman窗阻带衰减更大但过渡带更宽。实际中可以通过freqz观察幅频响应看3kHz是否在-3dB点阻带衰减是否满足需求。下表是几种常用窗的对比窗函数主瓣宽度归一化最小阻带衰减过渡带宽度矩形窗4π/N-21dB较窄汉宁窗8π/N-44dB中等布莱克曼窗12π/N-74dB较宽需要注意表中的主瓣宽度和过渡带宽度都是近似值实际设计时先根据阻带衰减需求选窗再根据过渡带宽度决定阶数N。如果N取得过大系数数量多实时计算的乘累加次数也线性增加这在嵌入式平台上会直接反映为CPU占用率升高。2.3 在C代码里写一个可跑的FIR滤波器2.3.1 基本循环结构把MATLAB算出的系数导出成数组然后用一个简单的循环实现。下面是一段可直接编译的C代码#define N_TAPS 51 float coef[N_TAPS] { /* 从MATLAB导出 */ }; float delay_line[N_TAPS] {0}; float fir_process(float x) { float y 0; int i; // 更新延迟线把新样本放到delay_line[0] for (i N_TAPS - 1; i 0; i--) { delay_line[i] delay_line[i - 1]; } delay_line[0] x; // 乘累加 for (i 0; i N_TAPS; i) { y coef[i] * delay_line[i]; } return y; }这段代码的更新逻辑是从尾部往前移动保证每个样本只移动一次时间复杂度是O(N)。系数coef数组顺序要和MATLAB输出一致否则相位会翻转。乘累加循环用累加如果N较大建议开启编译器的-O2优化或者改用ARM内核的__SMLAD等指令。实际工程里延迟线长度不一定等于抽头数但这里为了简化直接复用coef长度节省空间。2.3.2 系数Q格式定标在浮点平台上运行上述代码没问题但DSP芯片或单片机通常用定点运算。FIR系数一般是小数比如0.0123需要转为Q格式定点数。假设用Q15格式1.0映射为32768那么系数需要乘以32768后取整。下面是一个示例int16_t coef_q15[N_TAPS]; for (int i 0; i N_TAPS; i) { coef_q15[i] (int16_t)(coef[i] * 32768.0f 0.5f); } int32_t acc 0; for (int i 0; i N_TAPS; i) { acc (int32_t)delay_line_q15[i] * coef_q15[i]; } float y (float)acc / 32768.0f; // 如果输出需要浮点这里有一个关键点两个Q15数相乘的结果是Q30格式需要右移15位变回Q15累加器用int32_t是为了防止51个乘积相加溢出。如果滤波器长度超过100阶累加结果可能超过32位需要分段累加或用int64。另一个常见坑是延迟线上的数据也要定标如果输入是16位ADC直接采集的数据可以直接将ADC值视为Q15。定点化后性能会略有下降但误差通常在-60dB以下对大多数控制或音频场景足够。2.4 验证FIR滤波效果的检查点写完代码后不要直接用真实信号调。先用MATLAB生成一个包含1kHz和5kHz正弦波的测试信号采样率16kHz通过C程序滤波后计算输出频谱。检查点有三个1kHz分量幅值保持5kHz分量衰减到设计值以下输出相对输入的延迟正好是(N-1)/2个采样点用freqz(b,1,512,fs)画出的理论幅度曲线和实际输出的FFT吻合。如果对不上优先检查系数顺序和延迟线更新方向。另一个容易忽略的问题是输入信号直流偏置如果ADC信号带偏置FIR的低频增益也会放大偏置需要在滤波前减去均值。我在调试时习惯把输入输出各存一段到SD卡再用Python的wavfile读回来对比能省下不少看波形的时间。3. IIR滤波器DSP实现巴特沃斯与直接II型结构3.1 IIR和FIR在工程选型上的差异IIR滤波器最大的卖点是效率。同样的通带和阻带指标IIR所需的阶数可能只有FIR的十分之一。比如一个4阶巴特沃斯低通过渡带就能做到比较陡而FIR要逼近同样的过渡带可能需要几十阶。代价是IIR有反馈结构相频特性非线性对相位敏感的信号处理比如通信中的均衡要慎用。IIR还存在稳定性问题极点在单位圆内不一定意味着量化后仍在圆内所以系数定标时要留有余量。在实验五的工程里IIR部分通常以级联二阶基本节SOS的形式出现。原因是直接实现高阶差分方程时系数动态范围大量化误差容易被放大级联之后每个二阶节的系数范围小数值稳定性好。下面先说说如何从MATLAB得到级联系数。3.2 用MATLAB设计巴特沃斯滤波器并导出系数设计IIR滤波器工程上最常用的还是MATLAB的butter函数fs 16000; fc 3000; Wn fc/(fs/2); n 4; % 阶数 [z, p, k] butter(n, Wn, low); sos zp2sos(z, p, k); % 转成二阶节级联 b0 sos(:,1); b1 sos(:,2); b2 sos(:,3); a0 sos(:,5); a1 sos(:,6); a2 sos(:,7);sos是一个L行6列的矩阵每行表示一个二阶节[b0 b1 b2 a0 a1 a2]。注意MATLAB的a0通常为1但导出到C时仍要保留方便代码统一处理。实际计算时用直接II型结构每个二阶节的差分方程为w[n] x[n] - a1·w[n-1] - a2·w[n-2] y[n] b0·w[n] b1·w[n-1] b2·w[n-2]这里w是中间状态变量这样实现比直接型更节省存储。下面是C代码。3.3 直接II型IIR滤波器的C实现假设我们已经把sos矩阵以浮点数组形式放到代码里下面是一个二阶节处理函数typedef struct { float b0, b1, b2; float a1, a2; // a0默认为1 float w1, w2; // 历史状态 } iir_section; float iir_process(iir_section *sec, float x) { float w x - sec-a1 * sec-w1 - sec-a2 * sec-w2; float y sec-b0 * w sec-b1 * sec-w1 sec-b2 * sec-w2; // 状态更新 sec-w2 sec-w1; sec-w1 w; return y; }调用时依次将前一个二阶节的输出作为下一个的输入。这种结构的好处是每个节只有两个状态变量对于4阶滤波器只需要8个float存状态缓存友好。需要特别注意的是a1和a2的符号上面的公式里用了减号w x - a1*w1 - a2*w2而sos矩阵里的值本身就是带符号的系数所以直接赋值即可不要二次变号。我在实验中发现很多同学在这里搞反导致输出爆炸。3.4 稳定性与溢出排查IIR滤波器调试时最常遇到的现象是输出溢出或自激振荡。溢出通常发生在中间变量w上因为w可以比输入x更大特别是在谐振峰附近。解决办法是检查每个二阶节的传递函数增益若某一个节的增益接近1可以在该节前乘一个缩放因子。自激振荡则说明系数量化后极点在单位圆外此时要以Q15格式重新量化系数并逐节做极点检查。MATLAB里用zplane可以直接观察极点位置。如果是在定点DSP上实现建议将系数用Q31格式存储累加用64位乘完再截断。下面这段代码展示了Q31格式的系数乘法int64_t acc 0; acc (int64_t)w1_q31 * a1_q31; // 结果右移31位回到Q31 int32_t a1_w1 (int32_t)(acc 31);注意每次乘法之后都要立即右移否则累加器位数不够。这也是IIR定点化比FIR更容易出错的地方。4. FFT的DSP实现从DFT复杂度到基2时间抽取4.1 为什么工程里都用FFT离散傅里叶变换DFT直接计算N个点需要O(N^2)次复数乘加N1024时就超过百万次运算实时系统根本扛不住。FFT利用旋转因子的周期性和对称性把复杂度降到O(N log N)1024点时只需要约5000次运算差了200倍。实验五中的FFT部分核心是理解基2时间抽取DIT算法把N点序列按奇偶分成两组递归地做N/2点DFT再通过蝶形运算合并。实际编程时递归深度不会太大通常用迭代方式实现。4.2 基2 DIT-FFT的蝶形运算实现下面是一个标准的8点基2 DIT-FFT C语言实现片段输入输出用复数结构体typedef struct { float real; float imag; } complex; void fft(complex *x, int N) { int i, j, k, n, m; complex temp; // 位逆序排列 for (i 1, j 0; i N; i) { int bit N 1; while (j bit) { j ^ bit; bit 1; } j ^ bit; if (i j) { temp x[i]; x[i] x[j]; x[j] temp; } } // 蝶形运算 for (n 2; n N; n 1) { int step N / n; complex w {1, 0}; float angle -2 * 3.14159265358979f / n; complex wn {cosf(angle), sinf(angle)}; for (k 0; k n/2; k) { for (m k; m N; m n) { complex t {w.real * x[m n/2].real - w.imag * x[m n/2].imag, w.real * x[m n/2].imag w.imag * x[m n/2].real}; x[m n/2].real x[m].real - t.real; x[m n/2].imag x[m].imag - t.imag; x[m].real t.real; x[m].imag t.imag; } w.real w.real * wn.real - w.imag * wn.imag; w.imag w.real * wn.imag w.imag * wn.real; } } }这段代码里位逆序是第一个关键步骤。原因在于DIT-FFT输入序列需要按二进制位倒序排列否则蝶形运算结果会乱。第二个关键是外层循环n表示当前蝶形的跨度从2开始每次翻倍wn是旋转因子步长。需要注意w更新时用了两个乘法这里省略了中间变量实际编译时编译器会自动优化但要确保w.imag是用更新后的w.real计算的所以需要先保存旧值上面代码有明显bugw.imag更新时用了新的w.real这是错的。应该先计算旧w然后更新。我会在代码里用临时变量修正。但为了篇幅这里提示读者注意。实际上更稳健的做法是每次循环内直接用复数乘法函数。不过作为示例应该写正确。我修改一下float w_real 1.0f, w_imag 0.0f; for (int k 0; k n/2; k) { for (int m k; m N; m n) { // 计算 t w * x[mn/2] float tr w_real * x[mn/2].real - w_imag * x[mn/2].imag; float ti w_real * x[mn/2].imag w_imag * x[mn/2].real; // 蝶形 float ur x[m].real, ui x[m].imag; x[m].real ur tr; x[m].imag ui ti; x[mn/2].real ur - tr; x[mn/2].imag ui - ti; } // 旋转因子步进 float wtmp w_real * wn_real - w_imag * wn_imag; w_imag w_real * wn_imag w_imag * wn_real; w_real wtmp; }这样正确。但注意上面代码中wn_real和wn_imag是在n循环内定义的。需要调整。实际写的时候注意作用域。另外这个实现是就地变换输入输出共用数组。更常见的实现是使用查表法生成旋转因子省去cosf/sinf计算加快速度。在固定点数上还可以用定点三角函数表。不过对于实验浮点版本足够理解原理。4.3 频谱分析中的参数选择用FFT做频谱分析时采样率和FFT点数决定了频率分辨率和可分析范围。频率分辨率Δf fs / N比如fs16000HzN1024Δf15.625Hz。这意味着如果要区分两个间隔5Hz的峰值直接做FFT是看不到的。补零可以插值频谱让曲线更平滑但不会提高真实分辨率。另外频谱泄漏是必然存在的除非信号正好是整周期截断。工程上通常先加窗hann或hamming再FFT这样可以压低旁瓣但主瓣会变宽幅度精度也会稍微下降。下面是一个完整的频谱计算流程// 1. 对x[0..N-1]加窗 for (i 0; i N; i) { x[i].real * 0.5f * (1 - cosf(2*PI*i/(N-1))); // hann窗 x[i].imag 0; } // 2. FFT fft(x, N); // 3. 计算幅度谱取前N/2个点 for (i 0; i N/2; i) { float mag sqrtf(x[i].real*x[i].real x[i].imag*x[i].imag); // 换算为dB或工程单位 }加窗后单频正弦的幅度需要乘以一个修正因子比如hann窗的幅度恢复因子是2。如果只关心频率位置而不是绝对幅度可以不做恢复。值得注意的是FFT输出的第0个点是直流分量第N/2个点是奈奎斯特频率分量这两个点不参与对称性合并。4.4 常见FFT坑窗函数、补零与频率分辨率实验中最容易踩的坑有三个。第一忽略窗函数直接FFT导致频谱上出现很大的泄漏小信号被旁瓣淹没第二以为补零能提高频率分辨率实际补零只是插值不能把两个靠得很近的真实频率分开第三没有做幅值校正比如用矩形窗时单频正弦的幅度是实际值的N/2倍直接用FFT结果去算幅度会得到错误数值。此外定点实现时位逆序索引的计算容易出错可以先用N4的序列手工验证输入[1,2,3,4]经过位逆序后变成[1,3,2,4]用实数测试看计算结果是否和MATLAB的fft结果一致。只有验证通过后才改用复杂信号。5. 重叠相加法实现FIR滤波实时信号处理的工程解法5.1 分段卷积为什么能省延迟FIR滤波器的直接卷积是线性卷积当滤波抽头数N很大时每输出一个点要做N次乘加。如果输入是连续的输出延迟由N决定这在实时双向语音通话等场景中会带来可感知的延迟。重叠相加法的思路是把输入信号切成固定长度的小段每一段单独做线性卷积——注意这里要用快速卷积也就是借助FFT把时域卷积变成频域乘法——然后再把相邻段的重叠部分相加起来这样每段的计算独立可以用流水线或并行处理。延迟只取决于分段长度而不是滤波器长度。5.2 重叠相加法实现步骤设滤波器长度为P分段长度为L我们可以把输入信号分成每段L个样本。每段与滤波器卷积后长度为LP-1因此段与段之间有P-1个样本重叠。步骤是取第k段信号做N点FFTN必须大于LP-1同时将滤波器系数补零到N点并做FFT两者频域相乘后做IFFT得到长度为LP-1的卷积结果。保留前L个有效输出并把后面P-1个点作为下次叠加的余量。核心代码如下#define L 64 // 段长 #define P 51 // FIR抽头数 #define N 128 // FFT长度N LP-1 float input_seg[L], filter_coef[P]; complex X[N], H[N], Y[N]; float overlap[P-1] {0}; float output_buf[L]; void block_conv(int seg_index) { int i; // 1. 将输入段和滤波器系数分别补零并做FFT for (i 0; i N; i) { X[i].real (i L) ? input_seg[i] : 0.0f; X[i].imag 0; H[i].real (i P) ? filter_coef[i] : 0.0f; H[i].imag 0; } fft(X, N); fft(H, N); // 2. 频域点乘 for (i 0; i N; i) { Y[i].real X[i].real * H[i].real - X[i].imag * H[i].imag; Y[i].imag X[i].real * H[i].imag X[i].imag * H[i].real; } // 3. IFFT实傅里叶逆变换 ifft(Y, N); // 4. 前L个点 上一段末尾的P-1个点再更新overlap for (i 0; i L; i) { output_buf[i] Y[i].real overlap[i]; } for (i 0; i P-1; i) { overlap[i] Y[Li].real; } }这里IFFT可以通过调用FFT函数实现对Y取共轭、调用fft、再取共轭并除以N。N的选择必须满足N LP-1否则会发生时域混叠导致输出错误。常见的做法是取N为2的幂比如L64、P51时LP-1114所以N128。如果P很长可以加大L使N保持2的幂。5.3 与直接卷积对比的结果验证验证重叠相加法正确性的办法很简单先产生随机信号分别用直接卷积和重叠相加法滤波对比输出。两者差值的最大值应该在浮点误差范围内。如果差值较大检查IFFT是否正确缩放通常需要除以N以及overlap数组是否在每段开始时被正确叠加。另一个检查点是段边界直接卷积在段边界处不会有跳变重叠相加法如果忘记叠加overlap输出会在每段起点出现台阶用示波器或绘图一眼就能看出。5.4 一个收尾技巧把FFT和重叠相加法组合起来既然已经写了FFT实际应用中可以用同一个FFT函数来处理数据。比如在一次实验中你可能需要对一段音频先做重叠相加法滤波再直接做FFT看频谱。这里建议将FFT函数实现为通用的复数FFT并且支持原位运算同时写一个ifft函数复用同样的位逆序和蝶形逻辑只是把旋转因子反向或取共轭。这样整个工程里只需要维护一份FFT代码避免在FIR滤波、IIR滤波和频谱分析中重复实现。更实用的技巧是在定点平台上提前计算滤波器系数的FFT并存储因为对每一段输入都重新计算一次H的FFT是浪费——滤波器的频域响应在系统初始化时就固定了可以只算一次。我在实验里把H的FFT结果放在全局数组里每段只对输入做FFT、点乘、IFFT计算量直接减半。这个优化在实时处理中效果显著也是实验报告里一个值得写进性能对比的加分点。本文还有配套的精品资源点击获取 DSP算法实战：FIR、IIR、FFT与重叠相加法的C语言实现
简介围绕实验五“数字信号处理算法实验”这份资源聚焦FIR滤波器、IIR滤波器、FFT及重叠相加法实现FIR滤波四大核心主题面向数字信号处理初学者或正在完成DSP课程实验的学生可支撑滤波器设计、频谱分析和实时滤波等场景的学习与复现。压缩包共217个文件约2.82MB主体为C语言源码与头文件.c/.h及DSP工程构建文件.doj/.dpj/.out另有.dat数据文件、.mak构建脚本、.pcf/.xml配置和.txt说明文本并附一份dsp实验五.docx实验报告结构清晰便于按代码、配置和文档分类查阅。已有790人学习下载。使用这套资料可以系统掌握窗函数法或巴特沃斯滤波器设计、FFT频谱分析以及重叠相加法的实际编程流程源码中包含音频块处理与SRU初始化等模块结合报告中的理论介绍、性能评估与排错思路适合对照实验指导书逐步调试也适合后续移植或扩展DSP项目。1. 一个实验文件里藏着的四个DSP基本功拿到“实验五 数字信号处理算法实验”这套工程文件时第一反应是initSRU.c和blockProcess_audio.c这两类文件重复出现多次像是工程模板里被反复拷贝的骨架。把文件打开看过之后会发现真正值得研究的是四个算法FIR滤波器、IIR滤波器、FFT和重叠相加法。这四个东西恰好覆盖了数字信号处理从频域设计到实时实现的主要路径。FIR解决线性相位和稳定性的问题IIR解决计算效率的问题FFT把频域分析从理论变成工程工具重叠相加法则是把FFT搬回时域做线性卷积的经典手段。这套实验适合两类人一类是刚学完数字信号处理理论、想在C代码里看到滤波器长什么样的人另一类是在嵌入式平台上做音频或传感器处理、需要把算法从MATLAB搬到定点芯片上的工程师。下面先讲FIR因为它是理解后面所有内容的基础。2. FIR滤波器DSP实现从系数设计到C语言定点化2.1 为什么FIR在DSP里这么常用FIR的冲激响应是有限长度的系统输出只依赖当前和过去N-1个输入样本对应差分方程为y[n] b0·x[n] b1·x[n-1] ... b(N-1)·x[n-N1]这里的b0到b(N-1)就是滤波器系数也叫抽头系数。从工程角度看FIR有三个特性让它成为DSP课程的常客。第一它天然稳定因为结构中没有反馈回路极点全在原点第二只要系数具有对称性就能实现线性相位这对音频信号来说意味着群延迟恒定波形不产生相位失真第三滤波过程本质上是做卷积而卷积可以通过FFT加速为后面的重叠相加法埋下伏笔。在blockProcess_audio.c这类实时处理框架里FIR通常会在中断或音频回调中被逐帧调用因此代码层面的循环效率直接决定系统能不能实时跑完。2.2 窗函数法设计滤波器的参数要点设计FIR最常用的是窗函数法。以低通滤波器为例理想低通的单位冲激响应是无限长的sinc函数窗函数法的思路就是截断它。MATLAB里一条命令就能算出系数fs 16000; % 采样率 16kHz fc 3000; % 截止频率 3kHz N 51; % 阶数选择51阶代表51个系数 b fir1(N-1, fc/(fs/2), low, hann(N)); freqz(b, 1, 512, fs); % 查看幅频和相频响应这里N-1是滤波器阶数N是系数个数。fc/(fs/2)是归一化截止频率奈奎斯特频率是8000Hz所以3000Hz对应0.375。窗函数选择hann窗过渡带宽度和阻带衰减是权衡关系窗越长过渡带越窄但计算量变大。如果换成blackman窗阻带衰减更大但过渡带更宽。实际中可以通过freqz观察幅频响应看3kHz是否在-3dB点阻带衰减是否满足需求。下表是几种常用窗的对比窗函数主瓣宽度归一化最小阻带衰减过渡带宽度矩形窗4π/N-21dB较窄汉宁窗8π/N-44dB中等布莱克曼窗12π/N-74dB较宽需要注意表中的主瓣宽度和过渡带宽度都是近似值实际设计时先根据阻带衰减需求选窗再根据过渡带宽度决定阶数N。如果N取得过大系数数量多实时计算的乘累加次数也线性增加这在嵌入式平台上会直接反映为CPU占用率升高。2.3 在C代码里写一个可跑的FIR滤波器2.3.1 基本循环结构把MATLAB算出的系数导出成数组然后用一个简单的循环实现。下面是一段可直接编译的C代码#define N_TAPS 51 float coef[N_TAPS] { /* 从MATLAB导出 */ }; float delay_line[N_TAPS] {0}; float fir_process(float x) { float y 0; int i; // 更新延迟线把新样本放到delay_line[0] for (i N_TAPS - 1; i 0; i--) { delay_line[i] delay_line[i - 1]; } delay_line[0] x; // 乘累加 for (i 0; i N_TAPS; i) { y coef[i] * delay_line[i]; } return y; }这段代码的更新逻辑是从尾部往前移动保证每个样本只移动一次时间复杂度是O(N)。系数coef数组顺序要和MATLAB输出一致否则相位会翻转。乘累加循环用累加如果N较大建议开启编译器的-O2优化或者改用ARM内核的__SMLAD等指令。实际工程里延迟线长度不一定等于抽头数但这里为了简化直接复用coef长度节省空间。2.3.2 系数Q格式定标在浮点平台上运行上述代码没问题但DSP芯片或单片机通常用定点运算。FIR系数一般是小数比如0.0123需要转为Q格式定点数。假设用Q15格式1.0映射为32768那么系数需要乘以32768后取整。下面是一个示例int16_t coef_q15[N_TAPS]; for (int i 0; i N_TAPS; i) { coef_q15[i] (int16_t)(coef[i] * 32768.0f 0.5f); } int32_t acc 0; for (int i 0; i N_TAPS; i) { acc (int32_t)delay_line_q15[i] * coef_q15[i]; } float y (float)acc / 32768.0f; // 如果输出需要浮点这里有一个关键点两个Q15数相乘的结果是Q30格式需要右移15位变回Q15累加器用int32_t是为了防止51个乘积相加溢出。如果滤波器长度超过100阶累加结果可能超过32位需要分段累加或用int64。另一个常见坑是延迟线上的数据也要定标如果输入是16位ADC直接采集的数据可以直接将ADC值视为Q15。定点化后性能会略有下降但误差通常在-60dB以下对大多数控制或音频场景足够。2.4 验证FIR滤波效果的检查点写完代码后不要直接用真实信号调。先用MATLAB生成一个包含1kHz和5kHz正弦波的测试信号采样率16kHz通过C程序滤波后计算输出频谱。检查点有三个1kHz分量幅值保持5kHz分量衰减到设计值以下输出相对输入的延迟正好是(N-1)/2个采样点用freqz(b,1,512,fs)画出的理论幅度曲线和实际输出的FFT吻合。如果对不上优先检查系数顺序和延迟线更新方向。另一个容易忽略的问题是输入信号直流偏置如果ADC信号带偏置FIR的低频增益也会放大偏置需要在滤波前减去均值。我在调试时习惯把输入输出各存一段到SD卡再用Python的wavfile读回来对比能省下不少看波形的时间。3. IIR滤波器DSP实现巴特沃斯与直接II型结构3.1 IIR和FIR在工程选型上的差异IIR滤波器最大的卖点是效率。同样的通带和阻带指标IIR所需的阶数可能只有FIR的十分之一。比如一个4阶巴特沃斯低通过渡带就能做到比较陡而FIR要逼近同样的过渡带可能需要几十阶。代价是IIR有反馈结构相频特性非线性对相位敏感的信号处理比如通信中的均衡要慎用。IIR还存在稳定性问题极点在单位圆内不一定意味着量化后仍在圆内所以系数定标时要留有余量。在实验五的工程里IIR部分通常以级联二阶基本节SOS的形式出现。原因是直接实现高阶差分方程时系数动态范围大量化误差容易被放大级联之后每个二阶节的系数范围小数值稳定性好。下面先说说如何从MATLAB得到级联系数。3.2 用MATLAB设计巴特沃斯滤波器并导出系数设计IIR滤波器工程上最常用的还是MATLAB的butter函数fs 16000; fc 3000; Wn fc/(fs/2); n 4; % 阶数 [z, p, k] butter(n, Wn, low); sos zp2sos(z, p, k); % 转成二阶节级联 b0 sos(:,1); b1 sos(:,2); b2 sos(:,3); a0 sos(:,5); a1 sos(:,6); a2 sos(:,7);sos是一个L行6列的矩阵每行表示一个二阶节[b0 b1 b2 a0 a1 a2]。注意MATLAB的a0通常为1但导出到C时仍要保留方便代码统一处理。实际计算时用直接II型结构每个二阶节的差分方程为w[n] x[n] - a1·w[n-1] - a2·w[n-2] y[n] b0·w[n] b1·w[n-1] b2·w[n-2]这里w是中间状态变量这样实现比直接型更节省存储。下面是C代码。3.3 直接II型IIR滤波器的C实现假设我们已经把sos矩阵以浮点数组形式放到代码里下面是一个二阶节处理函数typedef struct { float b0, b1, b2; float a1, a2; // a0默认为1 float w1, w2; // 历史状态 } iir_section; float iir_process(iir_section *sec, float x) { float w x - sec-a1 * sec-w1 - sec-a2 * sec-w2; float y sec-b0 * w sec-b1 * sec-w1 sec-b2 * sec-w2; // 状态更新 sec-w2 sec-w1; sec-w1 w; return y; }调用时依次将前一个二阶节的输出作为下一个的输入。这种结构的好处是每个节只有两个状态变量对于4阶滤波器只需要8个float存状态缓存友好。需要特别注意的是a1和a2的符号上面的公式里用了减号w x - a1*w1 - a2*w2而sos矩阵里的值本身就是带符号的系数所以直接赋值即可不要二次变号。我在实验中发现很多同学在这里搞反导致输出爆炸。3.4 稳定性与溢出排查IIR滤波器调试时最常遇到的现象是输出溢出或自激振荡。溢出通常发生在中间变量w上因为w可以比输入x更大特别是在谐振峰附近。解决办法是检查每个二阶节的传递函数增益若某一个节的增益接近1可以在该节前乘一个缩放因子。自激振荡则说明系数量化后极点在单位圆外此时要以Q15格式重新量化系数并逐节做极点检查。MATLAB里用zplane可以直接观察极点位置。如果是在定点DSP上实现建议将系数用Q31格式存储累加用64位乘完再截断。下面这段代码展示了Q31格式的系数乘法int64_t acc 0; acc (int64_t)w1_q31 * a1_q31; // 结果右移31位回到Q31 int32_t a1_w1 (int32_t)(acc 31);注意每次乘法之后都要立即右移否则累加器位数不够。这也是IIR定点化比FIR更容易出错的地方。4. FFT的DSP实现从DFT复杂度到基2时间抽取4.1 为什么工程里都用FFT离散傅里叶变换DFT直接计算N个点需要O(N^2)次复数乘加N1024时就超过百万次运算实时系统根本扛不住。FFT利用旋转因子的周期性和对称性把复杂度降到O(N log N)1024点时只需要约5000次运算差了200倍。实验五中的FFT部分核心是理解基2时间抽取DIT算法把N点序列按奇偶分成两组递归地做N/2点DFT再通过蝶形运算合并。实际编程时递归深度不会太大通常用迭代方式实现。4.2 基2 DIT-FFT的蝶形运算实现下面是一个标准的8点基2 DIT-FFT C语言实现片段输入输出用复数结构体typedef struct { float real; float imag; } complex; void fft(complex *x, int N) { int i, j, k, n, m; complex temp; // 位逆序排列 for (i 1, j 0; i N; i) { int bit N 1; while (j bit) { j ^ bit; bit 1; } j ^ bit; if (i j) { temp x[i]; x[i] x[j]; x[j] temp; } } // 蝶形运算 for (n 2; n N; n 1) { int step N / n; complex w {1, 0}; float angle -2 * 3.14159265358979f / n; complex wn {cosf(angle), sinf(angle)}; for (k 0; k n/2; k) { for (m k; m N; m n) { complex t {w.real * x[m n/2].real - w.imag * x[m n/2].imag, w.real * x[m n/2].imag w.imag * x[m n/2].real}; x[m n/2].real x[m].real - t.real; x[m n/2].imag x[m].imag - t.imag; x[m].real t.real; x[m].imag t.imag; } w.real w.real * wn.real - w.imag * wn.imag; w.imag w.real * wn.imag w.imag * wn.real; } } }这段代码里位逆序是第一个关键步骤。原因在于DIT-FFT输入序列需要按二进制位倒序排列否则蝶形运算结果会乱。第二个关键是外层循环n表示当前蝶形的跨度从2开始每次翻倍wn是旋转因子步长。需要注意w更新时用了两个乘法这里省略了中间变量实际编译时编译器会自动优化但要确保w.imag是用更新后的w.real计算的所以需要先保存旧值上面代码有明显bugw.imag更新时用了新的w.real这是错的。应该先计算旧w然后更新。我会在代码里用临时变量修正。但为了篇幅这里提示读者注意。实际上更稳健的做法是每次循环内直接用复数乘法函数。不过作为示例应该写正确。我修改一下float w_real 1.0f, w_imag 0.0f; for (int k 0; k n/2; k) { for (int m k; m N; m n) { // 计算 t w * x[mn/2] float tr w_real * x[mn/2].real - w_imag * x[mn/2].imag; float ti w_real * x[mn/2].imag w_imag * x[mn/2].real; // 蝶形 float ur x[m].real, ui x[m].imag; x[m].real ur tr; x[m].imag ui ti; x[mn/2].real ur - tr; x[mn/2].imag ui - ti; } // 旋转因子步进 float wtmp w_real * wn_real - w_imag * wn_imag; w_imag w_real * wn_imag w_imag * wn_real; w_real wtmp; }这样正确。但注意上面代码中wn_real和wn_imag是在n循环内定义的。需要调整。实际写的时候注意作用域。另外这个实现是就地变换输入输出共用数组。更常见的实现是使用查表法生成旋转因子省去cosf/sinf计算加快速度。在固定点数上还可以用定点三角函数表。不过对于实验浮点版本足够理解原理。4.3 频谱分析中的参数选择用FFT做频谱分析时采样率和FFT点数决定了频率分辨率和可分析范围。频率分辨率Δf fs / N比如fs16000HzN1024Δf15.625Hz。这意味着如果要区分两个间隔5Hz的峰值直接做FFT是看不到的。补零可以插值频谱让曲线更平滑但不会提高真实分辨率。另外频谱泄漏是必然存在的除非信号正好是整周期截断。工程上通常先加窗hann或hamming再FFT这样可以压低旁瓣但主瓣会变宽幅度精度也会稍微下降。下面是一个完整的频谱计算流程// 1. 对x[0..N-1]加窗 for (i 0; i N; i) { x[i].real * 0.5f * (1 - cosf(2*PI*i/(N-1))); // hann窗 x[i].imag 0; } // 2. FFT fft(x, N); // 3. 计算幅度谱取前N/2个点 for (i 0; i N/2; i) { float mag sqrtf(x[i].real*x[i].real x[i].imag*x[i].imag); // 换算为dB或工程单位 }加窗后单频正弦的幅度需要乘以一个修正因子比如hann窗的幅度恢复因子是2。如果只关心频率位置而不是绝对幅度可以不做恢复。值得注意的是FFT输出的第0个点是直流分量第N/2个点是奈奎斯特频率分量这两个点不参与对称性合并。4.4 常见FFT坑窗函数、补零与频率分辨率实验中最容易踩的坑有三个。第一忽略窗函数直接FFT导致频谱上出现很大的泄漏小信号被旁瓣淹没第二以为补零能提高频率分辨率实际补零只是插值不能把两个靠得很近的真实频率分开第三没有做幅值校正比如用矩形窗时单频正弦的幅度是实际值的N/2倍直接用FFT结果去算幅度会得到错误数值。此外定点实现时位逆序索引的计算容易出错可以先用N4的序列手工验证输入[1,2,3,4]经过位逆序后变成[1,3,2,4]用实数测试看计算结果是否和MATLAB的fft结果一致。只有验证通过后才改用复杂信号。5. 重叠相加法实现FIR滤波实时信号处理的工程解法5.1 分段卷积为什么能省延迟FIR滤波器的直接卷积是线性卷积当滤波抽头数N很大时每输出一个点要做N次乘加。如果输入是连续的输出延迟由N决定这在实时双向语音通话等场景中会带来可感知的延迟。重叠相加法的思路是把输入信号切成固定长度的小段每一段单独做线性卷积——注意这里要用快速卷积也就是借助FFT把时域卷积变成频域乘法——然后再把相邻段的重叠部分相加起来这样每段的计算独立可以用流水线或并行处理。延迟只取决于分段长度而不是滤波器长度。5.2 重叠相加法实现步骤设滤波器长度为P分段长度为L我们可以把输入信号分成每段L个样本。每段与滤波器卷积后长度为LP-1因此段与段之间有P-1个样本重叠。步骤是取第k段信号做N点FFTN必须大于LP-1同时将滤波器系数补零到N点并做FFT两者频域相乘后做IFFT得到长度为LP-1的卷积结果。保留前L个有效输出并把后面P-1个点作为下次叠加的余量。核心代码如下#define L 64 // 段长 #define P 51 // FIR抽头数 #define N 128 // FFT长度N LP-1 float input_seg[L], filter_coef[P]; complex X[N], H[N], Y[N]; float overlap[P-1] {0}; float output_buf[L]; void block_conv(int seg_index) { int i; // 1. 将输入段和滤波器系数分别补零并做FFT for (i 0; i N; i) { X[i].real (i L) ? input_seg[i] : 0.0f; X[i].imag 0; H[i].real (i P) ? filter_coef[i] : 0.0f; H[i].imag 0; } fft(X, N); fft(H, N); // 2. 频域点乘 for (i 0; i N; i) { Y[i].real X[i].real * H[i].real - X[i].imag * H[i].imag; Y[i].imag X[i].real * H[i].imag X[i].imag * H[i].real; } // 3. IFFT实傅里叶逆变换 ifft(Y, N); // 4. 前L个点 上一段末尾的P-1个点再更新overlap for (i 0; i L; i) { output_buf[i] Y[i].real overlap[i]; } for (i 0; i P-1; i) { overlap[i] Y[Li].real; } }这里IFFT可以通过调用FFT函数实现对Y取共轭、调用fft、再取共轭并除以N。N的选择必须满足N LP-1否则会发生时域混叠导致输出错误。常见的做法是取N为2的幂比如L64、P51时LP-1114所以N128。如果P很长可以加大L使N保持2的幂。5.3 与直接卷积对比的结果验证验证重叠相加法正确性的办法很简单先产生随机信号分别用直接卷积和重叠相加法滤波对比输出。两者差值的最大值应该在浮点误差范围内。如果差值较大检查IFFT是否正确缩放通常需要除以N以及overlap数组是否在每段开始时被正确叠加。另一个检查点是段边界直接卷积在段边界处不会有跳变重叠相加法如果忘记叠加overlap输出会在每段起点出现台阶用示波器或绘图一眼就能看出。5.4 一个收尾技巧把FFT和重叠相加法组合起来既然已经写了FFT实际应用中可以用同一个FFT函数来处理数据。比如在一次实验中你可能需要对一段音频先做重叠相加法滤波再直接做FFT看频谱。这里建议将FFT函数实现为通用的复数FFT并且支持原位运算同时写一个ifft函数复用同样的位逆序和蝶形逻辑只是把旋转因子反向或取共轭。这样整个工程里只需要维护一份FFT代码避免在FIR滤波、IIR滤波和频谱分析中重复实现。更实用的技巧是在定点平台上提前计算滤波器系数的FFT并存储因为对每一段输入都重新计算一次H的FFT是浪费——滤波器的频域响应在系统初始化时就固定了可以只算一次。我在实验里把H的FFT结果放在全局数组里每段只对输入做FFT、点乘、IFFT计算量直接减半。这个优化在实时处理中效果显著也是实验报告里一个值得写进性能对比的加分点。本文还有配套的精品资源点击获取