DSP Question: invfreqs.m

1 次查看（过去 30 天）
Matlab2010
Matlab2010 2012-5-8
评论： Steven Lord 2019-9-27
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!
  1 个评论
Matlab2010
Matlab2010 2012-5-9
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

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Matlab2010
Matlab2010 2012-5-9
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');
  2 个评论
Harry
Harry 2019-9-27
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
Steven Lord 2019-9-27
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.

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Get Started with DSP System Toolbox 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by