DSP Question: invfreqs.m

John Michell answered . 2024-12-21 00:58:46

invfreqz.m has some odd rules re calling it. 

%%Orginal Data
  N = 5000;
  data = cumsum(randn(N,1));
  t = 252;
  a = 2 / (t+1);
  b = repmat(1-a,1 ,N).^(1:N); 
  b = b ./ sum(b); 
  a = 1;

%%Plot the Filter on some example data
ma = filter(b, a, data);
figure;plot(data); hold all; plot(ma, 'r');  

%%Plot the Response
figure;freqz(b,1, N);       
[h,w] = freqz(b,1,N);

%%Using dflit
%   b = 1; a = -1; %Sanity Check. simple difference filter
    % Hd = dfilt.df1([b a],1);   % num/ denom == a/b
    % fvtool(Hd);

%%Find the impluse response
% If n>(N-1) then you get a random answer! 
%If less then you get a visually different approximation,
%albiet it one with the same frequency response when plotted using freqz.m
n = N-1;  
m = 0;  % I choose 0 here as I have 1 in my orignal filter ==> the output comes out as aNew = 1;
wt = ones(size(w));

[bNew,aNew] = invfreqz(h,w,n,m);%, wt, 1000);  
%[bNew,aNew] = invfreqs(h,w,n,m);
%sys = tf(bNew,aNew);

%%Plot the filter coeffcients
% These now match if you used n = N-1
x1 = [0: 1/(size(b,2) -1)    : 1];
x2 = [0: 1/(size(bNew,2) -1) : 1];
figure;plot(x1,b); hold all; plot(x2,bNew, 'r');

%these always match
figure; freqz(b);
figure; freqz(bNew);
%invfreqz expects the complex-valued frequency response, not just the
%magnitude. If you only have the magnitude then use fdesign.arbmag
% http://www.mathworks.co.uk/products/dsp-system/demos.html?file=/products/demos/shipping/dsp/arbmagdemo.html
M = 1000; 
F = w/max(w); 
A = abs(h); 
d = fdesign.arbmag('N,F,A',M,F,A); 
Hd = design(d,'freqsamp'); 
%fvtool(Hd,'MagnitudeDisplay','Zero-phase');

b3 = Hd.Numerator(501:end); %symmetric so just take half. 
x1 = [0: 1/(size(b,2) -1) : 1]; 
x3 = [0: 1/(size(b3,2) -1) : 1]; 
figure;plot(x1,b); hold all; plot(x2,bNew, 'r'); plot(x3,b3, 'g');

Not satisfied with the answer ?? ASK NOW