FIR 滤波器用于频率的多频带滤波

Question

我的目标是像here那样进行多频带过滤，但是是在Python中。 Matlab 代码总结如下。与 Matlab 函数 kaiserord 非常相似的是 Python 实现 scipy.signal.kaiserord。但是，scipy.signal.kaiserord 不接受与 Matlab 版本相同的参数。

能否请有人提供一个 python 实现或关于如何在 python 中转换下面的代码的建议？

Fs = 100;                                                   % Sampling Frequency (Hz)
fcuts = [3.5 4 7 7.5 9.5 10 13 13.5 29 30 35 36];           % Frequencies
mags = [0 1 0 1 0 1 0 ];                                    % Passbands & Stopbands
devs = [0.05 0.01 0.05 0.01 0.05 0.01 0.05];                % Tolerances
[n,Wn,beta,ftype] = kaiserord(fcuts,mags,devs,Fs);          % Kaiser Window FIR Specification
n = n + rem(n,2);
hh = fir1(n,Wn,ftype,kaiser(n+1,beta),'noscale');           % Filter Realisation
figure
freqz(hh,1,2^14,Fs)
set(subplot(2,1,1), 'XLim',[0 50]);                         % Zoom Frequency Axis
set(subplot(2,1,2), 'XLim',[0 50]); 

Answer 1

我想我找到了办法。基于 Octave 代码here。下面是Matlab翻译。

功能：

function [n, w, beta, ftype] = kaiserOrd(f, m, dev, fs)
    if (nargin < 2 || nargin > 4)
        error('nargin < 2 || nargin > 4');
    end

    % Default sampling rate parameter
    if nargin < 4
        fs = 2;
    end

    % Parameter checking
    if length(f) ~= 2 * length(m) - 2
        error('kaiserord must have one magnitude for each frequency band');
    end

    if any(m(1:length(m) - 2) ~= m(3:length(m)))
        error('kaiserord pass and stop bands must be strictly alternating');
    end

    if length(dev) ~= length(m) && length(dev) ~= 1
        error('kaiserord must have one deviation for each frequency band');
    end

    dev = min(dev);

    if dev <= 0
        error('kaiserord must have dev > 0');
    end

    % Use midpoints of the transition region for band edges
    w = (f(1:2:length(f)) + f(2:2:length(f)))/fs;

    % Determine ftype
    if length(w) == 1
        if m(1) > m(2)
            ftype = 'LOW';
        else
            ftype = 'HIGH';
        end
    elseif length(w) == 2
        if m(1) > m(2)
            ftype = 'STOP';
        else
            ftype = 'BANDPASS';
        end
    else
        if m(1) > m(2)
            ftype = 'DC-1';
        else
            ftype = 'DC-0';
        end
    end

    % Compute beta from dev
    A = -20 * log10(dev);
    if (A > 50)
        beta = 0.1102 * (A - 8.7);
    elseif (A >= 21)
        beta = 0.5842 * (A - 21) ^ 0.4 + 0.07886 * (A - 21);
    else
        beta = 0.0;
    end

    % Compute n from beta and dev
    dw = 2 * pi * min(f(2:2:length(f)) - f(1:2:length(f))) / fs;
    n = max(1, ceil((A - 8) / (2.285 * dw)));

    % If the last band is high, make sure the order of the filter is even
    if ((m(1) > m(2)) == (rem(length(w), 2) == 0)) && rem(n, 2) == 1
        n = n + 1;
    end
end

演示：

clear all
close all

Fs = 1500;

duration  = 5;                    % second
N         = Fs*duration;          % sample length
t = 0:(1/(N-1)):1;                        % Time vector (1 second)

% Define the parameters for the sine waves
frequencies = [10, 50, 100, 150, 300, 400, 500, 1500]; % Sine wave frequencies (Hz)
amplitudes = [0.5, 0.8, 0.3, 0.1, 0.9, 1, 0.3, 1]; % Amplitudes of sine waves

% Initialize the signal
signal = zeros(size(t));

% Create the signal by adding multiple sine waves
for i = 1:length(frequencies)
    if mod(i,2) == 0
        signal = signal + amplitudes(i) * sin(2 * pi * frequencies(i) * t);
    else
        signal = signal + amplitudes(i) * cos(2 * pi * frequencies(i) * t);
    end
end

for i = 1:4
    figure
    i
    if i == 1
        bands = [300, 400];
        mag = [1, 0];
        dev = [0.1, 0.1];
    elseif i == 2
        bands = [400, 500] ;
        mag = [0, 1];
        dev = [0.1, 0.1];
    elseif i == 3
        bands = [100, 150, 400, 500];
        mag = [0, 1, 0];
        dev = [0.1, 0.1, 0.1];
    elseif i == 4
        bands = 10 * [10, 13, 15, 20, 30, 33, 35, 40] ;
        mag = [1, 0, 1, 0, 1];
        dev = [0.05, 0.05, 0.05, 0.05, 0.05];
    end
    
    % Use build in Matlab function
    % [n, w, beta, ftype] = kaiserord(bands, mag, dev, Fs);

    % Use hand made function
    [n, w, beta, ftype] = kaiserord(bands, mag, dev, Fs);

    d = max(1, fix(n / 10));
    if mag(length(mag)) == 1 && rem(d, 2) == 1
        d = d + 1;
    end

    [h, f] = freqz(fir1(n, w, ftype, kaiser(n + 1, beta), 'noscale'), 1, [], Fs);
    h = abs(h);
    hm = freqz(fir1(n - d, w, ftype, kaiser(n - d + 1, beta), 'noscale'), 1, [], Fs);
    
    subplot(3, 1, 1);
    plot(f, abs(hm), 'r', f, abs(h), 'b');
    title(sprintf('i=%d - order %d (red) - and order %d (blue)', i,  n - d, n));
    set(get(gca, 'Title'), 'FontSize', 8);

    b = [0, bands, Fs/2];
    hold on;
    for j = 2:2:length(b)
        hi = mag(j/2) + dev(1);
        lo = max(mag(j/2) - dev(1), 0);
        plot([b(j - 1), b(j), b(j), b(j - 1), b(j - 1)], [hi, hi, lo, lo, hi]);
    end
    hold off;
    
    % Apply the filter to the signal (remaining code is the same)
    signal_filtered = filtfilt(h, 1, signal);

    % Plot the original and filtered signals

    subplot(3, 1, 2);
    plot(t, signal);
    title('Original Signal');
    subplot(3, 1, 3);
    plot(t, signal_filtered);
    title('Filtered Signal');

end