CODEWITHSUNDEEP
pythonnumpyfiltersignal-processingfft
badaudio
badaudio

Reputation: 3

Numpy FFT error - Comb filter with envelope

I am currently having some issues understanding some discrepancies between the frequency response function calculated through the Z-transform and numpy's FFT algorithm. It is of a simple echo represented by the impulse response:

h[n] = δ[n] + αδ[n-L]

Where α is an attenuation factor for the echo and L is the delay amount in samples. The corresponding transfer function is given by:

H(f) = ( e^(j2πfΔL) + α ) / e^(j2πfΔL)

Where Δ is the sampling period.

I seem to be getting different results using the same number of frequency bins when directly plotting the transfer function magnitude above and when using numpy's fft algorithm. In particular, the FFT magnitude seems to form an envelope around the overall spectrum - I believe that I should be getting a simple comb filter as that of the transfer function method: imgur

Could anyone clarify why this may be happening and whether I have potentially overlooked anything? Is this due to errors in the DFT algorithms employed?

Appreciate your time, cheers!

import matplotlib.pyplot as pyplt
import numpy as np

fs = 48000  # Sample rate
T = 1/fs    # Sample period
L = 3000    # Delay
a = 0.5     # Attenuation factor

# h[n] = dirac[n] + a * dirac[n-L]
h = np.zeros(L)
h[0] = 1
h[L-1] = a

# Transfer function H
freqs = np.arange(0, fs, fs/(2*L))
e = np.exp(1j*freqs*2*np.pi*L*T)
H = (e + a)/(e)

# Transfer function H via FFT - same # of bins
H_FFT = np.fft.fft(h, 2*L)

pyplt.figure()
# Correct comb filter
pyplt.plot(np.abs(H))
# Runing FFT gives a form of strange envelope error
pyplt.plot(np.abs(H_FFT))
pyplt.legend()

Upvotes: 0

Views: 1647

Answers (1)

Aguy
Aguy

Reputation: 8059

Your code is almost OK. What you need to change:

  1. As far as I understand there is no reason the length of Fourier transform vector will be proportional to L. It needs to be the size of your sampling frequency, i.e. fs.

  2. Your L is too high. The result oscillates too quickly. Try with lower L.

Here is a modified code to show how it works, plotted in two different figures for clarity:

import matplotlib.pyplot as plt
import numpy as np

fs = 48000  # Sample rate
T = 1/fs    # Sample period
L = 3    # Delay
a = 0.5     # Attenuation factor

# h[n] = dirac[n] + a * dirac[n-L]
h = np.zeros(fs)
h[0] = 1
h[L] = a

# Transfer function H
freqs = np.arange(0, fs)
e = np.exp(1j*freqs*2*np.pi*L*T)
H = (e + a)/(e)

# Transfer function H via FFT - same # of bins
H_FFT = np.fft.fft(h)

# Correct comb filter
plt.figure()
plt.plot(np.abs(H))
# Runing FFT gives a form of strange envelope error
plt.figure()
plt.plot(np.abs(H_FFT))

Upvotes: 0

Related Questions