Making a Minimum-Phase FIR Filter in Python

# -*- coding: utf-8 -*- """ Created on Fri May 27 14:07:46 2016 @author: Judd """ #simulate downsampling from 88.2kHz to 44.1khz import numpy as np from scipy import signal import matplotlib.pyplot as plt fs_input=88200 fs_output=44100 # define parameters for a kaiser-windowed, linear-phase lowpasss filter: numtaps=511 kaiserbeta=20.0 min_nyquist=min(fs_input,fs_output)/2 cutoff = (min_nyquist * 21 / 22) / fs_input # display normalized cutoff frequency and cutoff frequency in Hz print(cutoff) print(cutoff*fs_input) # create the lowpass filter: lpass = signal.firwin(numtaps,cutoff,window=('kaiser',kaiserbeta),pass_zero=True,scale=True,nyq=0.5) # get frequency response of the lowpass filter: freq, response = signal.freqz(lpass) ampl = np.abs(response) # plot the frequency response of the linear-phase Lowpass filter (in Blue): fig1 = plt.figure() #plt.xlabel('frequency') plt.ylabel('dB') #plt.xscale('log') plt.axis([20500,24050,-250,6]) plt.grid(True) ax1 = fig1.add_subplot(111) ax1.plot(fs_input*freq/(2*np.pi), 20*np.log10(ampl), 'b-') # freq in Hz # convert Linear Phase Filter to Minimum Phase: fftsize = 32768 # needs to be fairly generous for good results ! maxphaseFilter=np.real(np.fft.ifft(np.exp(signal.hilbert(np.real(np.log(np.fft.fft(lpass,fftsize))))))) minphaseFilter=maxphaseFilter[::-1] # reverse finalFilt=minphaseFilter[0:numtaps-1] # get frequency response of minimum-phase filter: freq2, response2 = signal.freqz(minphaseFilter) ampl2 = np.abs(response2) # plot frequency response of Minimum-phase Filter (in Green) ax1.plot(fs_input*freq2/(2*np.pi), 20*np.log10(ampl2), 'g-') # freq in Hz # plot the actual impulse response (minphase in Green, linphase in Blue): fig2 = plt.figure() plt.axis([0,511,-1,1]) ax2=fig2.add_subplot(111) ax2.plot(lpass,'b-') ax2.plot(minphaseFilter,'g-') plt.show()