BUG: Different results between numpy.fft.rfft and scipy.signal.freqz

[happyTonakai]

Describe your issue.

I'm using the following two functions to calculate the frequency response of an impulse response: numpy.fft.rfft and scipy.signal.freqz, but I get different results. Let's simplify the problem using a small number of FFT.

import numpy as np
import scipy.signal as signal
h = np.array([1,2,3])
H1 = np.fft.rfft(h, n=4)
w, H2 = signal.freqz(h, 1, worN=len(H1), include_nyquist=True)

NFFT = 4, rfft calculates DFT at frequencies 0, pi/2, pi which is exactly equal to w returned by freqz

>>> print(w)
[0.         1.57079633 3.14159265]

However H2 is different from H1

>>> print(H1)
[ 6.+0.j -2.-2.j  2.+0.j]
>>> print(H2)
[ 6. +0.j          0.5-4.33012702j -1.5+0.8660254j ]

The DFT at pi/2 should be real-valued for a real-valued signal, so rfft gives the correct result.

From the source code of freqz in line 438-455 of scipy/signal/_filter_design.py we can see why it gives the wrong result. It doesn't give the nyquist result at all even include_nyquist is set to True.

Reproducing Code Example

import numpy as np
import scipy.signal as signal
h = np.array([1,2,3])
H1 = np.fft.rfft(h, n=4)
w, H2 = signal.freqz(h, 1, worN=len(H1), include_nyquist=True)
print(H1)
print(H2)

Error message

No error occurs.

SciPy/NumPy/Python version information

1.9.3 1.23.4 sys.version_info(major=3, minor=9, micro=0, releaselevel='final', serial=0)

[WarrenWeckesser]

Thanks for reporting the issue, @happyTonakai. That's definitely a bug. It works in some cases (specifically the IIR case where a is not a scalar, or when worN is less than the length of b), but in examples such yours or the following, setting include_nyquist=True doesn't actually cause the Nyquist frequency to be included in the calculation of h, so the values returned in h do not correspond to the frequencies returned in w (except at w=0):

In [56]: freqz([1, 2, 3, 2], worN=8)
Out[56]: 
(array([0.        , 0.39269908, 0.78539816, 1.17809725, 1.57079633,
        1.96349541, 2.35619449, 2.74889357]),
 array([ 8.        +0.j        ,  5.73444627-4.73444627j,
         1.        -5.82842712j, -2.20371254-3.20371254j,
        -2.        +0.j        , -0.03892814+1.03892814j,
         1.        +0.17157288j,  0.50819441-0.49180559j]))

In [57]: freqz([1, 2, 3, 2], worN=8, include_nyquist=True)
Out[57]: 
(array([0.        , 0.44879895, 0.8975979 , 1.34639685, 1.7951958 ,
        2.24399475, 2.6927937 , 3.14159265]),
 array([ 8.        +0.j        ,  5.73444627-4.73444627j,
         1.        -5.82842712j, -2.20371254-3.20371254j,
        -2.        +0.j        , -0.03892814+1.03892814j,
         1.        +0.17157288j,  0.50819441-0.49180559j]))

The parameter include_nyquist was added in #11368; @AtsushiSakai and @larsoner, can you take a look at this?

[larsoner]

Yikes I think we just changed the frequency returned, not the actual values.

To avoid this problem, I think we should do two things:

1. Don't use np.linspace to construct the frequencies, use rfftfreq instead. It makes sure they're always at least consistent.
2. Fix the actual bug mentioned here by changing the number of FFT points and slicing afterward correctly

@AtsushiSakai are you up for giving it a shot?

[drestebon]

I discussed this in #15273, there is a solution proposed in one my posts there.

The frequency grid is wrong, when using rfft. To illustrate, two calls of to freqz give different results:

import numpy as np
import scipy.signal as sig
d = sig.unit_impulse(4,3)
d_ = np.pad(d, (0,126-len(d)))
w, D = sig.freqz(d,1,worN=64,include_nyquist=True)
w, Dp = sig.freqz(d,1,worN=w)
D_ = np.fft.rfft(d_)
w_ = np.fft.rfftfreq(len(d_)) * 2 * np.pi
print('error freqz[rfft] v/s freqz[poly] = ', np.sum(np.abs(D-Dp)**2))
print('error freqz[rfft] v/s np.fft.rfft = ', np.sum(np.abs(D-D_)**2))
print('error freqz[poly] v/s np.fft.rfft = ', np.sum(np.abs(D_-Dp)**2))

The output:

error freqz[rfft] v/s freqz[poly] =  0.4657960974715796
error freqz[rfft] v/s np.fft.rfft =  0.4657960974715789
error freqz[poly] v/s np.fft.rfft =  1.1300975888209927e-29

[larsoner]

@drestebon would you be up for opening a PR to fix the accuracy / grid problem? Ultimately we should make sure the performance is correct first, i.e., in a first PR with test cases showing errant behavior on main and fixes on your PR. Then we can improve performance (assuming it's not already part of the fixes) in a separate PR. WDYT?

[drestebon]

@larsoner I will try to squeeze a couple of hours out of next week (or the following) ;)