Skip to content
Snippets Groups Projects
Commit b5a48dc4 authored by Eric Kooistra's avatar Eric Kooistra
Browse files

Filter design with fcutoff gain 0.5 or sqrt(0.5).

parent dfb2022f
Branches
No related tags found
No related merge requests found
Pipeline #84073 passed
Stage: linting
Show whitespace changes
Inline Side-by-side
applications/lofar2/model/pfb_os/filter_design_cutoff_frequency.ipynb 0 → 100644
+ 276
− 0
View file @ b5a48dc4
%% Cell type:markdown id:c69a2eb8 tags:
# Filter design with specific cutoff frequency
Author: Eric Kooistra, Jun 2024
Purpose:
* Practise DSP [1].
* Design filter with specific gain at cutoff frequency, typically gain 0.5 or sqrt(0.5)
References:
1. dsp_study_erko, summary of DSP books
%% Cell type:code id:02689e50 tags:
``` python
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
```
%% Cell type:code id:65235f50 tags:
``` python
# Auto reload module when it is changed
%load_ext autoreload
%autoreload 2
# Add rtdsp module path to $PYTHONPATH
import os
import sys
module_path = os.path.abspath(os.path.join('../'))
if module_path not in sys.path:
sys.path.insert(0, module_path)
# Import rtdsp
from rtdsp.firfilter import design_fir_low_pass_filter
from rtdsp.windows import half_cosine_window
from rtdsp.fourier import dtft, estimate_gain_at_frequency, estimate_frequency_at_gain
from rtdsp.plotting import plot_power_spectrum, plot_magnitude_spectrum
```
%% Output
The autoreload extension is already loaded. To reload it, use:
%reload_ext autoreload
%% Cell type:code id:c49515de tags:
``` python
# Filterbank
Ntaps = 16 # number of taps per poly phase FIR filter
Ndft = 1024 # DFT size
method = 'firls'
method = 'remez'
```
%% Cell type:code id:c5c90a6b tags:
``` python
# Samples
fs = 200e6 # sample rate
ts = 1 / fs # sample period
```
%% Cell type:code id:6d3a14bc tags:
``` python
# Subbands
Nsub = Ndft // 2 # number of subbands in fs / 2
fsub = fs / Ndft # subband frequency
tsub = 1 / fsub # subband period
```
%% Cell type:markdown id:15bf6804 tags:
# FIR low pass filter
%% Cell type:code id:0a69b385 tags:
``` python
# Specifications
Ncoefs = Ndft * Ntaps # use odd length for integer group delay
cutoffBeta = 0.2
fcutoff = fsub / 2 # center frequency of transition band
fpass = (1 - cutoffBeta) * fcutoff # pass band frequency
fstop = (1 + cutoffBeta) * fcutoff # stop band frequency
# fcutoff workaround for firls
tBW = (fstop - fpass) / 2
cutoffBW = tBW / 10 # is default value in design_fir_low_pass_filter()
# Gain at fcutoff center frequency of the transition band
cutoffGain = 0.5
#cutoffGain = np.sqrt(0.5)
# weight pass band ripple versus stop band ripple, and of zero width interval in the transition band
rippleWeights = [1, 1, 1]
hPrototype = design_fir_low_pass_filter(method,
Ncoefs, fpass, fstop, fcutoff, cutoffGain, rippleWeights, fs=fs, cutoffBW=cutoffBW)
if method == 'remez':
# Window tails to remove equiripple
hPrototype = hPrototype * half_cosine_window(Ncoefs, roBeta=0.1)
# Plot impulse response
plt.figure(1)
plt.plot(hPrototype, 'g')
plt.title('Impulse response')
plt.grid(True)
#plt.xlim((0, 200))
#plt.ylim((-0.00001, 0.00001))
# Plot transfer function (use frequency in fsub units)
h, f, HF = dtft(hPrototype)
plt.figure(2)
fLim = (-3, 3)
#fLim = None
dbLim = (-140, 5)
plot_power_spectrum(f, HF, 'g', fs / fsub, fLim, dbLim)
plt.xlabel('Frequency [fsub]')
plt.figure(3)
fLim = (0, 1) # Zoom in at cutoffGain
#fLim = (fpass / fsub, fstop / fsub) # Zoom in at cutoffGain
#fLim = None
voltLim = None
plot_magnitude_spectrum(f, HF, 'g', fs / fsub, fLim, voltLim)
plt.xlabel('Frequency [fsub]')
# Log abs(HF) at fcutoff
fIndex, fValue, fGain = estimate_gain_at_frequency(f * fs, HF, fcutoff)
print('estimate_gain_at_frequency fcutoff = %e [fsub]:' % (fcutoff / fsub))
print('. fIndex = %d' % fIndex)
print('. fValue = %e [fsub]' % (fValue / fsub))
print('. fGain = %e' % fGain)
fIndex, fValue, fGain = estimate_frequency_at_gain(f * fs, HF, cutoffGain)
print('estimate_frequency_at_gain cutoffGain = %e:' % cutoffGain)
print('. fIndex = %d' % fIndex)
print('. fValue = %e [fsub]' % (fValue / fsub))
print('. fGain = %e' % fGain)
```
%% Output
hInterpolated.imag ~= 0
> design_fir_low_pass_filter()
. Method = remez
. Q = 16.000000
. fpass = 78125.000000
. fstop = 117187.500000
. lpBW = 195312.500000
. fcutoff = 97656.250000
. cutoffGain = 0.500000
. fNyquist = 100000000.000000
. fs = 200000000.000000
. Ncoefs = 16384
. DC sum = 1.000000
. Symmetrical coefs = True
estimate_gain_at_frequency fcutoff = 5.000000e-01 [fsub]:
. fIndex = 131200
. fValue = 5.000000e-01 [fsub]
. fGain = 5.018223e-01
estimate_frequency_at_gain cutoffGain = 5.000000e-01:
. fIndex = 131200
. fValue = 5.000000e-01 [fsub]
. fGain = 5.018223e-01
%% Cell type:markdown id:e3542672 tags:
%% Cell type:code id:0b979a1f tags:
``` python
```
applications/lofar2/model/pfb_os/filter_design_erfc.ipynb 0 → 100644
+ 111
− 0
View file @ b5a48dc4
%% Cell type:markdown id:ea8e78bf tags:
# Filter design based on complementary error function (erfc)
Author: Eric Kooistra, Jun 2024
Purpose:
* Practise DSP [1].
* Design filter based on the complementary error function erfc = 1 - erf, as used in [2]
Description:
The erfc() is point symmeric around (0, 1) and decreases from 2 at -inf to 0 at +inf.
The shape of the erfc() is used as the template for the transfer function of a filter.
The advantage of using the erfc(), compared to e.g. a raised cosine filter with roll off beta, is that the shape has no discontinuities.
References:
1. dsp_study_erko, summary of DSP books
2. near perfect reconstruction, W. Lubberhuizen
%% Cell type:code id:7bf5ef09 tags:
``` python
import numpy as np
import matplotlib.pyplot as plt
from scipy import special
```
%% Cell type:code id:01918bc7 tags:
``` python
K = 4
Ndft = 16
Ntaps = 8
Ncoef = Ndft * Ntaps
F = np.arange(Ncoef) / Ncoef
x = K * (Ndft * F - 0.5)
A = 0.5 * special.erfc(x)
plt.plot(A)
plt.xlim((0,Ntaps))
print('%e' % A[0])
print('%e' % A[20])
print('%e' % A[-1])
```
%% Output
9.976611e-01
5.612149e-30
0.000000e+00
%% Cell type:code id:4611e655 tags:
``` python
```
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment