DSP Question: invfreqs.m

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!

[Matlab2010]

I try something else:

% http://www.mathworks.co.uk/products/dsp-system/demos.html?file=/products/demos/shipping/dsp/arbmagdemo.html

N = 1000;

F = w/max(w);

A = abs(h);

d = fdesign.arbmag('N,F,A',N,F,A);

Hd = design(d,'freqsamp');

%fvtool(Hd,'MagnitudeDisplay','Zero-phase');

b2 = Hd.Numerator(501:end); %symmetric so just take half.

x1 = [0: 1/(size(b,2) -1) : 1];

x2 = [0: 1/(size(b2,2) -1) : 1];

figure;plot(x1,b); hold all; plot(x2,b2, 'r');

Again, this is not even close

Answer 1

[Matlab2010]

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');

[Harry]

Thank you - I was starting to pull my hair out over this. I agree, the way to call it is very odd and it's badly documented by Matlab's usually high standards. Thanks for posting your clear example.

[Steven Lord]

If you have feedback on a function's documentation page, please click on one of the stars in the "How useful was this information?" box at the bottom of the page. Once you've done that, you can enter free-form text explaining why you chose that rating. In that text, you can express your confusion about what was written in the documentation, suggest improvements, or point out places where you believe it to be in error. The documentation staff reviews those ratings (I don't know the exact frequency of that review) and I know the documentation writers with whom I work have incorporated that feedback when they revise the pages they own.