matlab - Octave或Matlab中的invfreqz策略

Question

我有一个简单的低通滤波器的幅度和相位数据的集合。我可以使用 Octave 中的 invfreqz 成功地将滤波器拟合到复杂的频率响应中。但是，如果我希望以更高的采样率（比如 4 倍）使用相同的数据安装相同的滤波器怎么办？但是，没有高达新的 4 倍奈奎斯特频率的幅度和频率数据可用。无法在较高频率上收集磁力和相位数据，因此必须通过其他方法或近似值来添加。

将填充物幅度和相位数据“侵入”到新的奈奎斯特频率的最简单方法是什么，以便 invfreqz 有最好的机会以新的采样率实现对收集的数据的良好拟合？

Answer 1

这是我不久前编写的一个函数，可用于改变低/带通滤波器的采样率（我从来没有费心将其扩展到高通）。一般来说，您应该选择滤波器响应几乎为零的点作为“截止”频率（我不敢说“截止”，因为它不是正常意义上的截止）；高于此频率，响应设置为零。

function [b, a, fNew, HNew, fOld, HOld, HCut] = ...
    shiftFilterRate(b, a, fs1, fs2, order, fCut)
% SHIFTFILTERRATE   Shift sampling rate of filter
%
%   [bNew, aNew] = shiftFilterRate(b, a, fs1, fs2, order, fCut) designs a
%   FIR filter of the given order to operate at sampling frequency fs2, in
%   Hz, whose frequency magnitude characteristics match those of the filter
%   whose coefficients are b, a, that operates at sampling rate fs1, in Hz.
%   The function will only try to match the filter's magnitude in the
%   region [0 fCut], in Hz.
%
%   [bNew, aNew, fNew, HNew, fOld, HOld, HCut] = shiftFilterRate(...)
%   returns additional parameters:
%       fNew:   The frequencies at which the designed filter's transfer
%               function was evaluated (with resolution of 1 Hz)
%       HNew:   The transfer function of the designed filter
%       fOld:   The frequencies at which the existing filter's transfer
%               function was evaluated (with a resolution of 1 Hz)
%       HOld:   The transfer function of the existing filter.
%       HCut:   The desired frequency response of the filter used as input
%               to the function fir2.
%
%   FIXME: Make this work for high pass filters.
%

if nargin < 5, fCut = inf; end
if nargin < 4, order = 50; end

%% Zero padding in frequency domain
res = 1; % Hz resolution
N1 = fs1 / res; % points at resolution before padding

% Original freq response
[H, f] = freqz(b, a, N1);
nanInd = find(isnan(H));
% Stabilise. NOTE: Will break if nanInd contains last elem or there are
% multiple NaNs in a row
H(nanInd) = H(nanInd + 1); 
f = f / pi;

% Normalise cutoff freq
fCutNorm = fCut / (fs1 / 2);

% Cut frequencies above fCut, we don't really need them and it makes the
% FIR filter nasty
HCut = H;
HCut(f > fCutNorm) = 0;

% Create new freq response
NNew = ceil(fs2 / fs1 * length(HCut));
fNew = linspace(0, 1, NNew)';
HNew = [HCut; zeros(NNew - N1, 1)];

% Design filter
b = fir2(order, fNew, abs(HNew));
a = 1;

HOld = H;
fOld = f;

if nargout > 3
    HNew = freqz(b, a, length(fNew));
end

请注意，在您的情况下，您可能希望将其设置为 -20 dB re 1，因为这似乎是您的滤波器提供的最大衰减。在那一点上：它看起来不是一个很好的过滤器......你有什么理由必须匹配这个响应？您可能最好只设计一个具有相同截止频率的 Butterworth（您肯定会在阻带中获得更多衰减和更多线性相位）。