Python: High Pass FIR Filter with Windowing

Question

I would like to create a basic High Pass FIR filter by Windowing inside Python.

My code is below and intentionally idiomatic - I know that you can (most likely) complete this with a single line of code in Python, but I'm learning. I used the basic sinc function with a rectangular window: my output works for signals that are additive (f1 + f2) but not multiplicative (f1 * f2), where f1 = 25kHz and f2 = 1MHz.

My questions are: did I misunderstand something fundamental or is my code incorrect? In general, I would like to extract only the high pass signal (f2 = 1MHz) and filter out everything else. I also included screenshots of what is generated for (f1 + f2) and (f1 * f2):

import numpy as np import matplotlib.pyplot as plt # create an array of 1024 points sampled at 40MHz # [each sample is 25ns apart] Fs = 40e6 T = 1/Fs t = np.arange(0,(1024*T),T) # create an ip signal sampled at Fs, using two frequencies F_low = 25e3 # 25kHz F_high = 1e6 # 1MHz ip = np.sin(2*np.pi*F_low*t) + np.sin(2*np.pi*F_high*t) #ip = np.sin(2*np.pi*F_low*t) * np.sin(2*np.pi*F_high*t) op = [0]*len(ip) # Define - # Fsample = 40MHz # Fcutoff = 900kHz, # this gives the normalised transition freq, Ft Fc = 0.9e6 Ft = Fc/Fs Length = 101 M = Length - 1 Weight = [] for n in range(0, Length): if( n != (M/2) ): Weight.append( -np.sin(2*np.pi*Ft*(n-(M/2))) / (np.pi*(n-(M/2))) ) else: Weight.append( 1-2*Ft ) for n in range(len(Weight), len(ip)): y = 0 for i in range(0, len(Weight)): y += Weight[i]*ip[ni] op[n] = y plt.subplot(311) plt.plot(Weight,'ro', linewidth=3) plt.xlabel( 'weight number' ) plt.ylabel( 'weight value' ) plt.grid() plt.subplot(312) plt.plot( ip,'r-', linewidth=2) plt.xlabel( 'sample length' ) plt.ylabel( 'ip value' ) plt.grid() plt.subplot(313) plt.plot( op,'k-', linewidth=2) plt.xlabel( 'sample length' ) plt.ylabel( 'op value' ) plt.grid() plt.show() 

Answer 1

You misunderstood something fundamental. A filter with a sinc window is designed to separate linearly combined frequencies; that is, frequencies combined by addition, not frequencies combined by multiplication. See Chapter 5 , A Guide for Scientists and Engineers, Digital Signal Processing for more information.

Code based on scipy.signal will provide similar results to your code:

 from pylab import * import scipy.signal as signal # create an array of 1024 points sampled at 40MHz # [each sample is 25ns apart] Fs = 40e6 nyq = Fs / 2 T = 1/Fs t = np.arange(0,(1024*T),T) # create an ip signal sampled at Fs, using two frequencies F_low = 25e3 # 25kHz F_high = 1e6 # 1MHz ip_1 = np.sin(2*np.pi*F_low*t) + np.sin(2*np.pi*F_high*t) ip_2 = np.sin(2*np.pi*F_low*t) * np.sin(2*np.pi*F_high*t) Fc = 0.9e6 Length = 101 # create a low pass digital filter a = signal.firwin(Length, cutoff = F_high / nyq, window="hann") # create a high pass filter via signal inversion a = -a a[Length/2] = a[Length/2] + 1 figure() plot(a, 'ro') # apply the high pass filter to the two input signals op_1 = signal.lfilter(a, 1, ip_1) op_2 = signal.lfilter(a, 1, ip_2) figure() plot(ip_1) figure() plot(op_1) figure() plot(ip_2) figure() plot(op_2) 

Impulse Response:

Linear combined input:

Filtered Output:

Nonlinear combined input:

Filtered Output: