help-octave
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: freqz problems


From: Julius Smith
Subject: Re: freqz problems
Date: Tue, 4 Oct 2005 11:29:59 -0700

I believe the following patch to freqz.m fixes this:

*** freqz.m    2005/05/15 00:47:23    1.1
--- freqz.m    2005/10/04 18:22:24
***************
*** 143,149 ****
      hb = fft (postpad (b, extent));
      hb = hb(1:n);
    else
!     hb = polyval (postpad (b, k), exp (j*w));
    endif
    if (length (a) == 1)
      ha = a;
--- 143,149 ----
      hb = fft (postpad (b, extent));
      hb = hb(1:n);
    else
!     hb = polyval (postpad (b, k), exp (j*w)) .* exp(-j*(k-1)*w);
    endif
    if (length (a) == 1)
      ha = a;
***************
*** 151,157 ****
      ha = fft (postpad (a, extent));
      ha = ha(1:n);
    else
!     ha = polyval (postpad (a, k), exp (j*w));
    endif
    h = hb ./ ha;
    w = Fs * w / (2*pi);
--- 151,157 ----
      ha = fft (postpad (a, extent));
      ha = ha(1:n);
    else
!     ha = polyval (postpad (a, k), exp (j*w)) .* exp(-j*(k-1)*w);
    endif
    h = hb ./ ha;
    w = Fs * w / (2*pi);

---------------------------------

The proposed fix is equivalent to the perhaps clearer
hb = polyval (fliplr(postpad (b, k)), exp (-j*w));
etc., but should be faster.

On 9/28/05, Raman Venkataramani <address@hidden> wrote:
Hello,

I was wondering if this is a known bug. I am running octave version 2.1.50.
In this example, I am computing the frequency response of an FIR filter "h".

octave:96>
octave:96> h=[1, 0.5]
h =

  1.00000  0.50000

octave:97> H=freqz(h, 1, 512, "whole");
octave:98> H2=freqz(h, 1, 2*pi*[0:511]/512);
octave:99> norm(H-H2)
ans = 27.713
octave:100> norm(abs(H)-abs(H2))
ans =  5.3545e-15
octave:101>

I would have expected that H = H2, but only their magnitudes are equal. It
turns out that H is the correct answer and H2 is the reponse of a shifter
version of h with h(2) being the sample at the origin and h(1) is the
sampele at time -1.

Thanks,
Raman



-------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:   http://www.octave.org
How to fund new projects:  http://www.octave.org/funding.html
Subscription information:   http://www.octave.org/archive.html
-------------------------------------------------------------



reply via email to

[Prev in Thread] Current Thread [Next in Thread]