Search code examples
pythonfilteringsignal-processinglowpass-filter

Why does my windowed-sinc function have non-linear phase?


Similar to many tutorials on the web, I've tried implementing a windowed-sinc lowpass filter using the following python functions:

def black_wind(w):
''' blackman window of width w'''
    samps = np.arange(w)
    return (0.42  - 0.5 * np.cos(2 * np.pi * samps/ (w-1)) + 0.08 * np.cos(4 * np.pi * samps/ (w-1)))

def lp_win_sinc(tw, fc, n):
''' lowpass sinc impulse response 
Parameters:
    tw = approximate transition width [fraction of nyquist freq]
    fc = cutoff freq [fraction of nyquest freq]
    n = length of output. 
Returns:
    s = impulse response of windowed-sinc filter appended zero-padding
    to make len(s) = n
'''
    m = int(np.ceil( 4./tw / 2) * 2) 
    samps = np.arange(m+1)
    shift = samps - m/2
    shift[m/2] = 1
    h = np.sin(2 * np.pi * fc * shift)/shift
    h[m/2] = 2 * np.pi * fc
    h = h * black_wind(m+1)
    h = h / h.sum()
    s = np.zeros(n)
    s[:len(h)] = h
    return s

For input: 'tw = 0.05', 'fc = 0.2', 'n = 6000', the magnitude of the fft seems reasonable.

tw = 0.05
fc = 0.2
n = 6000
lp = lp_win_sinc(tw, fc, n)
f_lp = np.fft.rfft(lp)
plt.figure()
x = np.linspace(0, 0.5, len(f_lp))
plt.plot(x, np.abs(f_lp))

magnitude of lowpass filter response

however, the phase is non-linear above ~fc.

plt.figure()
x = np.linspace(0, 0.5, len(f_lp))
plt.plot(x, np.unwrap(np.angle(f_lp)))

phase of lowpass filter response

Given the symmetry of the non-zero-padded portion of the impulse response, I would expect the resulting phase to be linear. Can someone explain what is going on? Perhaps I'm using a numpy function incorrectly, or maybe my expectations are incorrect. I'm very grateful for any help.

***********************EDIT***********************

based on some of the helpful comments to this question and some more work, I wrote a function that produces zero phase delay and is therefore a bit easier to interpret the np.angle() results.

def lp_win_sinc(tw, fc, n):
    m = int(np.ceil( 2./tw) * 2) 
    samps = np.arange(m+1)
    shift = samps - m/2
    shift[m/2] = 1
    h = np.sin(2 * np.pi * fc * shift)/shift
    h[m/2] = 2 * np.pi * fc
    h = h * np.blackman(m+1)
    h = h / h.sum()
    s = np.zeros(n)
    s[:len(h)] = h
    return np.roll(s, -m/2)

The main change here is using np.roll() to place the line of symmetry at t=0.


Solution

  • The magnitudes in the stop band are crossing zero. The phase of the coefficient after a zero crossing will jump by 180 degrees, which is confusing np.angle()/np.unwrap(). -1*180° = 1*0°