MATLAB: DSP Question: invfreqs.m

invfreqs

I generate some coefficents for a filter and can inspect the frequency response as following:
%%Orginal Data
N = 5000;
data = cumsum(randn(N,1));
t = 252;
a = 2 / (t+1);
b = repmat(1-a,1 ,N).^(1:N); %b are your filter coeff 
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);          
[h,w] = freqz(b,1);
I now explain my problem. I am now in the situation where I have a frequency response (i.e. the vector "h") and know nothing else.
I would like to estimate from this my original "b" (the filter coefficents) to allow me to estimate my variable "t".
I thought I could use invfreqs.m (or invfreqz.m) to do this, but Im afraid I dont know how.
%%Find the impluse response
n = 10;  % I choose a large number allowing a good approximation
m = 0;  % I choose 0 here as I have 1 in my orignal filter ==> the output comes out as aNew = 1;
[bNew,aNew] = invfreqz(h,w,n,m);
%[bNew,aNew] = invfreqs(h,w,n,m);
sys = tf(bNew,aNew)
%%Plot the filter coeffcients
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');
When I inspect the final plot, I would expect to see the red line (bNew) as a good approximation to b. It is not. not even close.
Clearly I am doing something very wrong. Please could someone with experince of how this function works, explain my mistake.
many thanks!

Best Answer

Solved. invfreqz.m has some odd rules re calling it. Could do with better documentation/examples to my mind.
%%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');
Related Question