Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Signal #6

Merged
merged 10 commits into from
Jun 12, 2024
Merged

Signal #6

merged 10 commits into from
Jun 12, 2024

Conversation

AngelEzquerra
Copy link
Contributor

This adds a signal processing related module, starting with a kaiser window and a firls FIR design function.

…nction)

This module will have signal processing related functions (e.g. to generate windows, resampling, etc).

The first set of functions added to this module allow you to calculate the Kaiser window.
To do so, we also add an implementation of the Bessel function of the first kind of order 0 (i0).

The code in this commit is based on the code for the corresponding functions in the scipy library.
…ign)

Adds a new `firls` procedure that can be used to design a linear phase FIR filter using least-squares error criterion. This function is avaiable in `scpy` (as `scipy.signal.firls`). This version is similar but has some differences in its interface and also as some additional functionality (it can design a broader set of filter types).

`firls` eturns the coefficients of the linear phase FIR filter of the requested
order and type which best fullfills the desired frequency response
"constraints" described by the `bands` (i.e. frequency), `desired` (i.e. gain)
and `weights` (e.i. constraint weights) input tensors.

Depending on the order and the value of the `symmetric` argument, the
output filter coefficients will correspond to a Type I, II, III or IV
FIR filter.

The constraints are specified as a pair of tensors describing "constrained
frequency bands" indicating the start and end frequencies of each band `k`,
as well as the gain of the filter at those edge frequencies.

Inputs:
  - fir_order: The "order" of the filter (which is 1 less than the length
               of the output tensor).
  - bands: 2-column matrix of non-overlapping, normalized frequency band
           edges in ascending order. Each row in the matrix corresponds to
           the edges of a constrained frequency band (the column contains
           the start frequency of the band and the second column contains
           the end frequency).
           Frequencies between the specified bands are "unconstrained"
           (i.e. ignored during the error minimization process).
           Frequency values must be in the [0.0, fs/2] range (i.e. they
           cannot be negative or exceed the Nyquist frequency).
  - desired: 2-column matrix that specifies the gains of the desired
             frequency response at the band "edges". Thus the lengths of
             the `bands` and `desired` tensors must match. If they do not,
             an exception is raised.
             For each band `k` the desired filter frequency
             response is such that its gain linearly changes from the
             `desired[k, 0]` value at the start of the band (i.e. at the
             `bands[k, 0]` frequency) to the value `desired[k, 1]` at the
             end of the band (i.e. at `bands[k, 1]`).
  - weights: Optional rank-1 Tensor of weights. Controls which frequency
             response "contraints" are given more "weight" during the
             least-squares error minimization process. The default is that
             all constraints are given the same weight. If provided, its
             length must be half the length of `bands` (i.e. there must be
             one constraint per band). An exception is raised otherwise.
  - symmetric: When `true` (the default), the result will be a symmetric
               FIR filter (Type I when `fir_order` is even and Type II
               when `fir_order` is odd). When `false`, the result will be
               an anti-symmetric FIR filter (Type III when `fir_order` is
               even and Type IV when `fir_order` is odd).
  - fs: The sampling frequency of the signal (as a float). Each frequency
        in `bands` must be between 0.0 and `fs/2` (inclusive).
        Default is 2.0, which means that by default the band frequencies
        are expected to be on the 0.0 to 1.0 range.

Result:
  - A Tensor containing the `fir_order + 1` taps of the FIR filter that
    best approximates the desired frequency response (i.e. the filter that
    minimizes the least-squares error vs the given constraints).

Notes:
  - Contrary to numpy's firls, the first argument is the FIR filter order,
    not the FIR length. The filter length is `filter_order + 1`.
  - Frequencies between "constrained bands" are considered "don't care"
    regions for which the error is not minimized.
  - When designing a filter with a gain other than zero at the Nyquist
    frequency (i.e. when `bands[^1] != 0.0`), such as high-pass and
    band-stop filters, the filter order must be even. An exception is
    raised otherwise.
  - The `bands` and `desired` can also be flat rank-1 tensors, as long as
    their length is even (so that they can be reshaped into 2 column
    matrices.

Examples:
```nim
echo firls(4, [[0.0, 0.3], [0.4, 1.0]].toTensor, [[1.0, 1.0], [0.0, 0.0]].toTensor)

echo firls(4, [0.0, 0.3, 0.4, 1.0].toTensor, [1.0, 1.0, 0.0, 0.0].toTensor)

echo firls(4,
  [[0.0, 1.5], [2.0, 5.0]].toTensor,
  [[1.0, 1.0], [0.0, 0.0]].toTensor,
  fs = 10.0)

echo firls(4,
  [[0.0, 0.3], [0.4, 1.0]].toTensor,
  [[1.0, 1.0], [0.0, 0.0]].toTensor,
  weights = [1.0, 0.5].toTensor)

echo firls(5,
  [[0.0, 0.3], [0.3, 0.6], [0.6, 1.0]].toTensor,
  [[1.0, 1.0], [1.0, 0.2], [0.0, 0.0]].toTensor)

echo firls(6,
  [[0.0, 0.4], [0.6, 1.0]].toTensor, [[0.0, 0.0], [0.9, 1.0]].toTensor,
  symmetric = false)

Trying to design a high-pass filter with odd order generates an exception:
echo firls(5,
  [0.0, 0.5, 0.6, 1.0].toTensor, [0.0, 0.0, 1.0, 1.0].toTensor,
  symmetric = false)
```
This procedure, which is available both in Matlab and scipy (as `scipy.signal.upfirdn`)  upsamples a rank-1 tensor, applies a FIR filter and downsamples it. It is basically a shortcut for a combination of upsample and convolve with a downsample factor. However this is a core signal processing operation that it deserves its own function.
Also update the description a little.
@AngelEzquerra AngelEzquerra mentioned this pull request Jun 4, 2024
@Vindaar
Copy link
Member

Vindaar commented Jun 7, 2024

For the Windows CI we need to install BLAS / LAPACK. Should be an easy extension of the CI YAML file. In ggplotnim I do it like this:

https://github.com/Vindaar/ggplotnim/blob/master/.github/workflows/ci.yml#L57-L71

We'll probably want to switch away from setup-nim soon, too. It hasn't been updated in a long time and I think at this point there are better ways to build / get Nim (for the OSX failure).

@Vindaar
Copy link
Member

Vindaar commented Jun 7, 2024

Otherwise looking good, as far as my code review can judge (i.e. haven't tried to check the implementation does what it says it does 😁). Only wondered if there's some refactoring that could be done to share more code between the symmetric / non symmetric branches in firls. There seems to be a lot of similarity in the for loop code. But I don't mind either way.

@AngelEzquerra
Copy link
Contributor Author

Otherwise looking good, as far as my code review can judge (i.e. haven't tried to check the implementation does what it says it does 😁). Only wondered if there's some refactoring that could be done to share more code between the symmetric / non symmetric branches in firls. There seems to be a lot of similarity in the for loop code. But I don't mind either way.

I'd be OK trying to unify the code a bit because as you said it both branches are quite similar. However, the differences are large enough that it might be a bit complicated, I'd have to give it a try.
If you don't mind, I'd rather add this as is, and then I can try to improve it later?

This adds a couple of `resample` procedures as well as some useful, supporting functions.
These procedurs let you resample a tensor by a certain up / down sampling ratio.
One version of the `resample` procedure lets you provide a specific antialiasing filter that will be used between the upsampling and downsampling operations. Another version automatically calculates an appropriate filter (using `firls`), given a specific "filter order factor" and a beta value for the kaiser window that is applied to the designed filter.
@Vindaar
Copy link
Member

Vindaar commented Jun 10, 2024

Will you add the BLAS / LAPACK dependcny for Windows into the CI file? Fixing OSX can be done at a later time. Once I look into it, I need to do it for a variety of repositories anyway.

@Vindaar
Copy link
Member

Vindaar commented Jun 12, 2024

Merging despite the CI issues. They are not very relevant imo for the time being.

@Vindaar Vindaar merged commit f7418dd into SciNim:master Jun 12, 2024
3 of 9 checks passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants