过滤多个频段

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

figure
Fs = 11025;

t = 0:1/Fs:1;              % Time vector (1 second)

% Define the parameters for the sine waves
frequencies = [4, 10, 30]; % Sine wave frequencies (Hz)
amplitudes = [0.5, 0.8, 0.3]; % Amplitudes of sine waves

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

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

for i = 1:4
    figure
    i
    if i == 1
        subplot(2, 2, 1);
        bands = [1200, 1500];
        mag = [1, 0];
        dev = [0.1, 0.1];
    elseif i == 2
        subplot(2, 2, 2);
        bands = [1000, 1500];
        mag = [0, 1];
        dev = [0.1, 0.1];
    elseif i == 3
        subplot(2, 2, 3);
        bands = [1000, 1200, 3000, 3500];
        mag = [0, 1, 0];
        dev = [0.1, 0.1, 0.1];
    elseif i == 4
        subplot(2, 2, 4);
        bands = 100 * [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('order %d (red) and order %d (blue)', 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

您可以使用 Matlab 中的构建函数或手工制作的函数 (kaiserOrd) 来运行演示。输出是相同的：