Add complex valued bandpass filter #468

zsoerenm · 2022-03-30T10:24:47Z

This pull request allows to create ComplexBandpass filters for one sided pass bands (different pass bands for positive or negative frequencies)
I haven't implemented scalefactor(coefs::Vector, ftype::ComplexBandpass) yet, as I'm not really sure how to implement it for this case. So this currently only works for scale = false.
Also I haven't implemented tests, yet. Where do all the filter coefficients in the tests for the different filters come from? Maybe someone can help me out with this.
Here is an example:

using DSP, FFTW, Plots
signal = randn(ComplexF64, 10000)
filtered_signal = filt(digitalfilter(ComplexBandpass(-5e6, -3e6; fs = 10e6), FIRWindow(hanning(512), scale = false)), signal)
p = periodogram(filtered_signal, fs = 10e6)
Plots.plot(fftshift(freq(p)), 10 * log10.(fftshift(power(p))))

martinholters · 2022-06-21T05:55:20Z

src/Filters/design.jl

@@ -241,6 +241,11 @@ function normalize_freq(w::Real, fs::Real)
    f >= 1 && error("frequencies must be less than the Nyquist frequency $(fs/2)")
    f
 end
+function normalize_complex_freq(w::Real, fs::Real)
+    f = 2*w/fs
+    f >= 2 && error("frequencies must be less than the Nyquist frequency $(fs)")


Should we also put a lower bound on f? It seems that 0 <= f < 2 and -1 <= f < 1 would be reasonable limits, but depending on context, one or the other may be more useful. So require -1 <= f < 2? Looks rather arbitrary. Or maybe we put no restriction on f at all? Really, I don't know.

martinholters · 2022-06-21T05:59:11Z

Regarding tests, I think most of the reference data has been produced using comparable Matlab functions. Personally, I'm not too huge a fan of this. For the case at hand, I could imagine determining the frequency responses of a couple of designed filters and verifying that they fulfill the specifications (under some assumptions on what passband ripple and stopband gain can be reasonably expected).

codecov · 2024-02-05T17:45:40Z

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 97.76%. Comparing base (76001f4) to head (1257934).

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #468      +/-   ##
==========================================
+ Coverage   97.74%   97.76%   +0.01%     
==========================================
  Files          18       18              
  Lines        3199     3225      +26     
==========================================
+ Hits         3127     3153      +26     
  Misses         72       72

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

martinholters

I think constructor should return values of the type they're asked to construct, which appears not to be the case here.

src/Filters/design.jl

wheeheee · 2024-04-30T22:16:58Z

I think this PR just needs a reference signal for the test, probably written with scipy (replacing the placeholder one I left). Not sure how the complex bandpass is supposed to work, so I'm unable to do that...

reorganize

martinholters · 2024-05-21T13:44:51Z

src/Filters/design.jl

-    w::T
+function normalize_complex_freq(w::Real, fs::Real)
+    f = 2 * w / fs
+    f >= 2 && throw(DomainError(f, "frequencies must be less than the Nyquist frequency $(fs)"))


Suggested change

f >= 2 && throw(DomainError(f, "frequencies must be less than the Nyquist frequency $(fs)"))

f >= 2 && throw(DomainError(f, "frequencies must be less than the sampling frequency $(fs)"))

martinholters · 2024-05-21T13:46:38Z

LGTM modulo the minor docstring fix. Good to go then or are there any objections?

wheeheee · 2024-05-21T14:11:27Z

scalefactor isn't implemented for complex bandpass. Is that meaningful or is it ok as is now?

martinholters · 2024-05-22T06:10:57Z

Ah, good point. I think this should work:

function scalefactor(coefs::Vector{T}, ftype::ComplexBandpass, fs::Real) where T<:Number
    n = length(coefs)
    freq = normalize_complex_freq(middle(ftype.w1, ftype.w2), fs)
    c = zero(T)
    for k = 1:n
        c = muladd(coefs[k], cispi(-freq * (k - (n + 1) / 2)), c)
    end
    return abs(c)
end

martinholters

Trying to re-create the example from the OP, I got different results which lead me to what I suspect to be a bug in how fs is used.

src/Filters/design.jl

martinholters · 2024-05-22T06:29:55Z

test/filter_design.jl

+    winfirtaps_jl = digitalfilter(ComplexBandpass(0.1, 0.2), FIRWindow(hamming(128), scale=false); fs=1)
+    winfirtaps_scipy = readdlm(joinpath(dirname(@__FILE__), "data", "digitalfilter_hamming_128_complexbandpass_fc1.txt"), '\t', ComplexF64)
+    @test winfirtaps_jl ≈ winfirtaps_scipy  # replace with real test data


Can we add a comment how the reference data was generated here? Also, I think we should check with another fs, just to make sure the scaling works as expected, e.g.

digitalfilter(ComplexBandpass(100, 200), FIRWindow(hamming(128), scale=false); fs=1000)

should give the same result.

I wrote # TODO: write actual tests for ComplexBandpass and # replace with real test data in the commit "add placeholder tests" to indicate dummy reference data, generated by the same function it was supposed to be tested against. Probably wasn't a good idea, but I intended to set it up so the reference data could just be changed to something legitimate afterwards...

Ah, I wasn't aware of this circulus vitiosus. (Even noting the comment, the _scipy suffix in the variable name lead me to wrong conclusions.)

I'll try to cook up a more effective test.

These tests differ from those of the other `digitalfilter` types in that they don't compare to a reference impulse response, but rather check the frequency response for plausibility.

martinholters reviewed Jun 21, 2022

View reviewed changes

wheeheee added this to the 0.8 milestone Feb 8, 2024

wheeheee force-pushed the ss/complex-bandpass-filter branch from 08c30a6 to 72a188a Compare February 14, 2024 05:51

wheeheee force-pushed the ss/complex-bandpass-filter branch from 9fd0547 to 5e30601 Compare February 26, 2024 10:51

wheeheee force-pushed the ss/complex-bandpass-filter branch from 5e30601 to 9a2d9e6 Compare March 8, 2024 15:57

martinholters requested changes Mar 11, 2024

View reviewed changes

src/Filters/design.jl Outdated Show resolved Hide resolved

wheeheee force-pushed the ss/complex-bandpass-filter branch from 9a2d9e6 to 70119bf Compare March 11, 2024 11:14

wheeheee reviewed Mar 11, 2024

View reviewed changes

src/Filters/design.jl Show resolved Hide resolved

martinholters reviewed Mar 11, 2024

View reviewed changes

src/Filters/design.jl Show resolved Hide resolved

wheeheee force-pushed the ss/complex-bandpass-filter branch from 70119bf to 4941d9d Compare April 26, 2024 21:26

zsoerenm and others added 5 commits May 14, 2024 22:30

Add complex valued bandpass filter

f97010b

reorganise constructors, cispi

cf5e4f5

reorganize

add placeholder tests

9cac7a1

DomainError

efd3c2d

fix constructors

373f26f

wheeheee force-pushed the ss/complex-bandpass-filter branch from b7cb383 to 373f26f Compare May 14, 2024 20:30

martinholters reviewed May 21, 2024

View reviewed changes

martinholters requested changes May 22, 2024

View reviewed changes

martinholters added 4 commits May 24, 2024 08:21

Add scalefactor for ComplexBandpass

0f8d2d8

Remove dummy test for ComplexBandpass

6023427

Fix frequency scaling in firprototype for ComplexBandpass

7a5efdb

Add tests for ComplexBandpass

1257934

These tests differ from those of the other `digitalfilter` types in that they don't compare to a reference impulse response, but rather check the frequency response for plausibility.

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Add complex valued bandpass filter #468

Add complex valued bandpass filter #468

zsoerenm commented Mar 30, 2022

martinholters Jun 21, 2022

martinholters commented Jun 21, 2022

codecov bot commented Feb 5, 2024 •

edited

Loading

martinholters left a comment

wheeheee commented Apr 30, 2024

martinholters May 21, 2024

martinholters commented May 21, 2024

wheeheee commented May 21, 2024

martinholters commented May 22, 2024

martinholters left a comment

martinholters May 22, 2024

wheeheee May 22, 2024

martinholters May 24, 2024

	f >= 2 && throw(DomainError(f, "frequencies must be less than the Nyquist frequency $(fs)"))
	f >= 2 && throw(DomainError(f, "frequencies must be less than the sampling frequency $(fs)"))

Add complex valued bandpass filter #468

Are you sure you want to change the base?

Add complex valued bandpass filter #468

Conversation

zsoerenm commented Mar 30, 2022

martinholters Jun 21, 2022

Choose a reason for hiding this comment

martinholters commented Jun 21, 2022

codecov bot commented Feb 5, 2024 • edited Loading

Codecov Report

martinholters left a comment

Choose a reason for hiding this comment

wheeheee commented Apr 30, 2024

martinholters May 21, 2024

Choose a reason for hiding this comment

martinholters commented May 21, 2024

wheeheee commented May 21, 2024

martinholters commented May 22, 2024

martinholters left a comment

Choose a reason for hiding this comment

martinholters May 22, 2024

Choose a reason for hiding this comment

wheeheee May 22, 2024

Choose a reason for hiding this comment

martinholters May 24, 2024

Choose a reason for hiding this comment

codecov bot commented Feb 5, 2024 •

edited

Loading