Python：窗口化高通FIR滤波器 [英] Python: High Pass FIR Filter by Windowing

问题描述

我的代码位于下方，故意惯用 - 我知道你可以（很可能）用Python中的一行代码来完成这个工作，但我正在学习。我已经使用了一个矩形窗口的基本sinc函数：我的输出适用于加法（f1 + f2），但不乘法（f1 * f2），其中f1 = 25kHz和f2 = 1MHz的信号。

我的问题是：我误解了一些基本的东西，或者是我的代码错了？
总之，我想只提取高通信号（f2 = 1MHz），并过滤其他的东西。我还包括了为（f1 + f2）和（f1 * f2）生成的内容截图： stack.imgur.com/DKtvx.pngalt =

$ p $ import numpy as np
import matplotlib.pyplot as plt

＃创建一个以40MHz采样的1024点阵列
＃[每个采样间隔25ns]
Fs = 40e6
T = 1 / F $
t = np.arange（0，（1024 * T），T）

＃创建一个以Fs采样的ip信号，使用两个频率
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，
＃这给出了归一化的转换频率，Ft
Fc = 0.9e6
Ft = Fc / Fs
长度= 101
M =长度 - 1
重量= []
在范围内（0，长度）：
if （n！=（M / 2））：
Weight.append（-np.sin（2 * np.pi * Ft *（n-（M / 2）））/（np.pi *（n-（M / 2））））
else： （len（Weight），len（ip））：
Weight.append（1-2 * Ft）





y = 0
在范围内（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（）


解决方案

你误解了一些基本的东西。加窗正弦滤波器被设计成分离线性组合频率;即通过加法组合的频率，而不是通过乘法组合的频率。有关详细信息，请参阅科学家和工程师指南的第5章至
数字信号处理。 。

基于scipy.signal的代码将提供与您的代码类似的结果：

 从pylab导入* 
导入scipy.signal作为信号
 
＃创建以40MHz采样的1024点阵列
＃[每个样本相距25ns] 
 Fs = 40e6 
 nyq = Fs / 2 
 T = 1 / Fs 
t = np.arange（0，（1024 * T），T）
 
＃创建一个以Fs采样的ip信号，使用两个频率：
 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 
长度= 101 
 
＃创建低通数字滤波器
a = signal.firwin（Length，cutoff = F_high / nyq， window =hann）
 
＃通过信号反转创建一个高通滤波器
a = -a 
a [Length / 2] = a [Length / 2] + 1 
 $ b $ figure（）
 plot（a，'ro'）
 
＃将高通滤波器应用于两个输入信号
 op_1 = signal.lfilter（一个， 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）
   
   
 
 pngalt =Impulse Response> 
 
 
 
线性组合输入：
 
 
  
 
 
 
已过滤的输出： 
 
 
  
 
 
 
非 - 非线性组合输入： 
 
 
  
 
 
 
已过滤的输出： 
 
 
  
 
I'd like to create a basic High Pass FIR Filter by Windowing within Python.

My code is below and is intentionally idiomatic - I'm aware you can (most likely) complete this with a single line of code in Python but I'm learning. I have used a basic a 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: Have I misunderstood something fundamental or is my code wrong? 
In summary, I'd like to extract just the high pass signal (f2=1MHz) and filter everything else out. I've also included screen shots 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[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()

解决方案
You've misunderstood something fundamental. The windowed sinc filter is designed to separate linearly combined frequencies; i.e. frequencies combined through addition, not frequencies combined through multiplication. See chapter 5 of The Scientist and Engineer's Guide to
Digital Signal Processing for more details. 

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:



Linearly Combined Input:



Filtered Output:



Non-linearly Combined Input:



Filtered Output:



                        
这篇关于Python：窗口化高通FIR滤波器的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！