Python：通過開窗高通FIR濾波器

Question

我想在Python中通過開窗創建一個基本的高通FIR濾波器。

我的代碼在下面，故意慣用 - 我知道你可以（很可能）用Python中的一行代碼來完成它，但我正在學習。我用矩形窗口使用了一個基本的sinc函數：我的輸出適用於加法（f1 + f2）但不乘法（f1 * f2）的信號，其中f1 = 25kHz，f2 = 1MHz。

我的問題是：我誤解了一些基本的或我的代碼錯誤？總之，我只想提取高通信號（f2 = 1MHz）並過濾其他所有信號。我還包括爲（F1 + F2）和（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[n-i] 
    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

你誤會了一些基本的東西。開窗式sinc濾波器設計用於分離線性組合頻率;即通過加法組合的頻率，而不是通過乘法組合的頻率。有關更多詳細信息，請參見「科學家和工程師指南 數字信號處理」的chapter 5。基於scipy.signal

代碼將提供類似的結果代碼：

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) 

衝激響應：

線性組合輸入：

過濾ED輸出：

非線性組合輸入：

濾波輸出：