CMSIS FIR 带通滤波器答案

【问题标题】：CMSIS FIR bandpass filterCMSIS FIR 带通滤波器
【发布时间】：2016-07-08 10:28:21
【问题描述】：

我正在尝试在 STM32F407 微控制器上实现 60kHz 带通滤波器，但遇到了一些问题。我在 MATLAB fdatool 的帮助下生成了过滤器，然后也在 MATLAB 中对其进行了仿真。以下 MATLAB 脚本对其进行了模拟。

% FIR Window Bandpass filter designed using the FIR1 function.
% All frequency values are in Hz.
Fs = 5250000;    % Sampling Frequency

N    = 1800;     % Order
Fc1  = 59950;    % First Cutoff Frequency
Fc2  = 60050;    % Second Cutoff Frequency
flag = 'scale';  % Sampling Flag
% Create the window vector for the design algorithm.
win = hamming(N+1);

% Calculate the coefficients using the FIR1 function.
b  = fir1(N, [Fc1 Fc2]/(Fs/2), 'bandpass', win, flag);
Hd = dfilt.dffir(b);
%----------------------------------------------------------
%----------------------------------------------------------
T = 1 / Fs;          % sample time
L = 4500;            % Length of signal
t = (0:L-1)*T;       % Time vector

% Animate the passband frequency span
for f=55500:50:63500
    signal = sin(2*pi*f*t);
    plot(filter(Hd, signal));
    axis([0 L -1 1]);

    str=sprintf('Signal frequency (Hz) %d', f);
    title(str);
    drawnow;
end

pause;
close all;

signal = sin(2*pi*50000*t) + sin(2*pi*60000*t) + sin(2*pi*78000*t);
signal = signal / 3;
signal = signal(1:1:4500);

filterInput = signal;
filterOutput = filter(Hd,signal);

subplot(2,1,1);
plot(filterInput);
axis([0 4500 -1 1]);

subplot(2,1,2);
plot(filterOutput)
axis([0 4500 -1 1]);
pause;

close all;

我从 fdatool 中将滤波器系数提取为 q15 格式的 16 位无符号整数，这是因为我使用的是 12 位 ADC。 MATLAB 生成的滤波器系数标头为here，得到的系数图如下图

下面是过滤器实现的代码，它显然不起作用，我真的不知道我能做些什么不同，我在网上查看了一些示例 Example 1 和 Example 2

#include "fdacoefs.h"

#define FILTER_SAMPLES 4500
#define BLOCK_SIZE     900    

static uint16_t firInput[FILTER_SAMPLES];
static uint16_t firOutput[FILTER_SAMPLES];
static uint16_t firState[NUM_TAPS + BLOCK_SIZE - 1];

uint16_t util_calculate_filter(uint16_t *buffer, uint32_t len)
{
    uint16_t i;   
    uint16_t max;
    uint16_t min;
    uint32_t index;

    // Create filter instance
    arm_fir_instance_q15 instance; 

    // Ensure that the buffer length isn't longer than the sample size
    if (len > FILTER_SAMPLES)
        len = FILTER_SAMPLES;   

    for (i = 0; i < len ; i++) 
    {
        firInput[i] = buffer[i];        
    }

    // Call Initialization function for the filter 
    arm_fir_init_q15(&instance, NUM_TAPS, &firCoeffs, &firState, BLOCK_SIZE);

    // Call the FIR process function, num of blocks to process = (FILTER_SAMPLES / BLOCK_SIZE)
    for (i = 0; i < (FILTER_SAMPLES / BLOCK_SIZE); i++) // 
    {
        // BLOCK_SIZE = samples to process per call
        arm_fir_q15(&instance, &firInput[i * BLOCK_SIZE], &firOutput[i * BLOCK_SIZE], BLOCK_SIZE); 
    }

    arm_max_q15(&firOutput, len, &max, &index);    
    arm_min_q15(&firOutput, len, &min, &index);

    // Convert output back to uint16 for plotting
    for (i = 0; i < (len); i++) 
    {
        buffer[i] = (uint16_t)(firOutput[i] - 30967);
    }

    return (uint16_t)((max+min));
}

ADC 以 5.25 MSPS 采样，它对 60kHz 信号进行 4500 次采样，在这里您可以看到滤波器的 Input，然后是滤波器的 Output，这很奇怪..

有什么明显的我错过了吗？因为我完全迷路了，任何指示和提示都是有帮助的！

【问题讨论】：

在 32 位 MCU 上进行如此多的 uint16_t 计算，我会立即怀疑各种溢出和隐式提升错误。您的代码有几个地方依赖于隐式整数提升。这通常表明程序员不考虑/不了解隐式提升，这反过来会产生非常微妙的错误。
为什么在 32 位 MCU 上使用 uint16_t 计算很重要？有处理半字的指令，q15 值也表示为 16 位无符号值，ADC 将 12 位样本存储在无符号 16 位缓冲区中，那么类型转换发生在哪里？
32 位 MCU 意味着您的 int 类型是 32 位，这反过来意味着您使用的所有小整数类型将被隐式提升为 32 位有符号类型。这是 C 语言强制要求的，与您的特定处理器或编译器无关。您需要阅读整数提升规则和平衡（“通常的算术转换”）如何工作。 Here's an explanation from CERT-C.
我知道，但我不明白这可能是什么问题？
变量在计算前改变符号。变量比较给出了意想不到的结果。看似有效但隐藏的潜在值溢出的计算会在轻微的代码更改时突然暴露。按位运算给出非常奇怪的结果。等等。很多潜在的细微错误。

标签： c matlab filter arm signal-processing

【解决方案1】：

正如 Lundin 指出的那样，我将其更改为使用 32 位整数，这实际上解决了我的问题。当然，我使用 MATLABS fdatool 作为有符号 32 位整数生成了新的滤波器系数。

static signed int firInput[FILTER_SAMPLES];
static signed int firOutput[FILTER_SAMPLES];
static signed int firState[NUM_TAPS + BLOCK_SIZE -1];

uint16_t util_calculate_filter(uint16_t *buffer, uint32_t len)
{
    uint16_t i;   
    int power;
    uint32_t index;

    // Create filter instance
    arm_fir_instance_q31 instance; 

    // Ensure that the buffer length isn't longer than the sample size
    if (len > FILTER_SAMPLES)
        len = FILTER_SAMPLES;   

   for (i = 0; i < len ; i++) 
    {
        firInput[i] = (int)buffer[i];        
    }

    // Call Initialization function for the filter 
    arm_fir_init_q31(&instance, NUM_TAPS, &firCoeffs, &firState, BLOCK_SIZE);

    // Call the FIR process function, num of blocks to process = (FILTER_SAMPLES / BLOCK_SIZE)
    for (i = 0; i < (FILTER_SAMPLES / BLOCK_SIZE); i++) // 
    {
        // BLOCK_SIZE = samples to process per call
        //arm_fir_q31(&instance, &firInput[i * BLOCK_SIZE], &firOutput[i * BLOCK_SIZE], BLOCK_SIZE); 
        arm_fir_q31(&instance, &firInput[i * BLOCK_SIZE], &firOutput[i * BLOCK_SIZE], BLOCK_SIZE); 
    }

    arm_power_q31(&firOutput, len, &power);

    // Convert output back to uint16 for plotting
    for (i = 0; i < (len); i++) 
    {
        buffer[i] = (uint16_t)(firOutput[i] - 63500);
    }

    return (uint16_t)((power/10));
}

【讨论】：