strategy for invfreqz in Octave or Matlab

Question

I have a collection of magnitude and phase data of a simple low pass filter. I can successfully fit a filter to the complex freq response using invfreqz in Octave. However, what if I wish to fit the same filter at a higher sampling rate ,say 4x, with the same data? However, there is no magnitude and freq data available up to the new 4x nyquist frequency. Mag and phase data cannot be collected at the higher frequencies so it must be added by some other method or approximation.

What would be the easiest way to "hack-in" filler magnitude and phase data up to the new nyquist frequency so that invfreqz would have the best possible chance of achieving a good fit for the collected data at the new sample rate?

What kind of filter is it? If it's low-pass or band-pass, this is easy. Most high-pass filters should be fine as well. Can you put an image of the magnitude / phase plot on tinypic? Essentially you want to pick some range up to which you're going to match the response, then design a new filter with fir2 or similar. Details will follow in an answer if you post a picture of the response. — wakjah
Hi. Thanks for response. It is lowpass filter. I have placed image of mag and phase on tinypic. tinypic.com/r/b7yxxl/6 and tinypic.com/r/n4e8pi/6 — user2243673
The data goes up to 7214. I can fit with invfreq with Fs=14428. However, what if I wish to use the data as a filter with much higher sample rate, say, 14428*4=57712? I need to fudge data up to 28,856. What is best way to do this? — user2243673

wakjah wakjah · Accepted Answer · 2013-04-05T09:16:25

Here is a function I wrote a while ago that can be used for shifting sampling rates of low/band-pass filters (I never bothered to extend it to high-pass). In general you should pick the point where your filter's response drops to almost nothing as the 'cut' frequency (I hesitate to say 'cut-off', because it is not a cut-off in the normal sense); above this frequency the response is set to zero.

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

Note that in your case, you might want to set it to -20 dB re 1, as this appears to be the maximum attenuation your filter provides. On that point: it doesn't look like it's a very good filter... Is there any reason why you have to match this response? You might be better off just designing, e.g., a Butterworth with the same cut-off (you'd certainly get more attenuation in the stop band and more linear phase).

strategy for invfreqz in Octave or Matlab

1 Answers