Skip to content

Transforms

Jump to bottom
Tim Sharii edited this page Jul 25, 2019 · 9 revisions

Currently, the following basic digital transforms are implemented:

Also read this:

Note, the implementation of algorithms is far from optimized. Nevertheless, even on humble modern hardware the overall speed of transforms allows performing a lot of useful things without significant or noticeable delays. And who says that it's the last version of NWaves? 😃 There's always a place for optimizations, and corresponding PRs are very welcome.

The FFT code is conventional decimation radix-2 FFT algorithm (O(Nlog(N))). Hence, direct/inverse FFT requires the length of input data to be a power of 2, otherwise exception will be thrown.

The FWT is coded according to classic Mallat decomposition / reconstruction algorithm.

Implementation of DCT I-IV transformers is straightforward (O(N^2)), since so far it's used only in MFCC-like algorithms where small array sizes are involved.

Other transforms are based on FFT, so their performance is conditioned by the performance of FFT.

For each transform there's a corresponding transformer object. Each transformer object has Direct() and Inverse() methods.

FFT transformers

There are two classes providing FFT functionality:

  • Fft - FFT for complex input
  • RealFft - FFT for real input

They give same results for real inputs. Most likely, you'll be working with real inputs, so in most practical cases RealFft should be preferred, since it performs significantly faster (about 25-35% gain). As of ver.0.9.2 all feature extractors, convolvers and filters use RealFft objects under the hood for computations.

The only parameter in Fft constructor is the size of FFT (that should be power of 2, as mentioned above). Both direct and inverse transforms are performed in-place on arrays of real and imaginary parts of samples. If the length of data differs from specified size of Fft object, exception will be thrown.

var fft = new Fft(1024);

float[] real = signal.First(1024).Samples;
float[] imag = new float [1024];

// in-place FFT
fft.Direct(real, imag);

// ...do something with real and imaginary parts of the spectrum...

// in-place IFFT
fft.Inverse(real, imag);

RealFft's analogous methods have slightly different prototypes. The following illustration explains the difference:

    Fft direct  : complex -> complex
RealFft direct  : real    -> complex
    Fft inverse : complex -> complex
RealFft inverse : complex -> real

So the code given above would change to:

var fft = new RealFft(1024);

// in-place FFT: real  ->   (real, imag)
fft.Direct(real, real, imag);

// ...do something with real and imaginary parts of the spectrum...

// in-place IFFT: (real, imag)  ->   real
fft.Inverse(real, imag, real);


// if original array must be preserved, another block can be used for output:

fft.Direct(input, outRe, outIm);
fft.Direct(outRe, outIm, output);

All other functions in Fft and RealFft classes are completely identical.

In some cases it's OK to work directly with real and imag arrays (see the code of Phase Vocoder and its variations). But often what we need is to compute some real-valued post-processed results of complex FFT, like magnitude spectrum or power spectrum:

var magnitudeSpectrum = 
    fft.MagnitudeSpectrum(signal[1000, 2024]);

var powerSpectrum = 
    fft.PowerSpectrum(signal.First(1024), normalize: false);

var logPowerSpectrum = 
    fft.PowerSpectrum(signal.Last(1024))
       .Samples
       .Select(s => Scale.ToDecibel(s))
       .ToArray();

The second parameter defines whether the spectrum should be normalized by the size of FFT or not. By default power spectrum is normalized, and magnitude spectrum is not normalized.

Also, methods MagnitudeSpectrum and PowerSpectrum will automatically zero-pad or truncate samples to fit FFT size (so unlike Direct() and Inverse() methods there will be no exceptions in this case).

The magnitude/power spectrum has length (FFT size) / 2 + 1, i.e. it includes 0th coefficient (DC component) and does not contain the symmetric part of the spectrum.

Finally, methods MagnitudeSpectrum and PowerSpectrum have overloaded versions for float[] arrays designed for efficient sequential spectral analysis of a signal (discussed below).

Fft.Shift is quite conventional routine in DSP. It shifts the zero-frequency component of spectrum to the center of the spectrum. It can work with any even-sized array:

float[] array = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15 };

Fft.Shift(array);

// array = { 8, 9, 10, 11, 12, 13, 14, 15, 0, 1, 2, 3, 4, 5, 6, 7 }

Short-Time Fourier Transform:

var stft = new Stft(1024, 250);
var timefreq = stft.Direct(signal);
var reconstructed = stft.Inverse(timefreq);

The first parameter in the constructor is the window size (actually, it may be less than FFT size and doesn't have to be a power of 2, however in ver.0.9.2 only windowSize = fftSize is allowed because of a little bug which will be fixed in newer versions; anyway, it's not crucial since usually there's no need to make windowSize not equal to fftSize in STFT). The second parameter is the hop size that can be any arbitrary positive value. The third parameter is the type of a window function. By default Hann window is used due to its "good" mathematical properties, so don't change it unless you do some exotic processing. The fourth parameter is the size of FFT that is calculated by default as the nearest power of 2 to window size.

The result of direct STFT has type List<Tuple<float[], float[]>>, i.e. collection of spectra. Each spectrum is represented with two arrays containing real and imaginary parts of spectral components. This collection contains all necessary data for inverse STFT. Often this information is redundant, and we need to evaluate only the so called spectrogram of a signal which is essentially a sequence of power spectra. In this case we'd use the Spectrogram method that returns simply the List<float[]>:

var stft = new Stft(512, 256);
var spectrogram = stft.Spectrogram(signal);

Sometimes it's convenient to work with STFT results in Magnitude-Phase form:

var stft = new Stft(1024, 250);
var magPhase = stft.MagnitudePhaseSpectrogram(signal);

// magnitudes and phases in the first frame:

float[] magnitudes = magPhase[0].Magnitudes;
float[] phases = magPhase[0].Phases;

// magPhase is of type MagnitudePhaseList:

//    public struct MagnitudePhaseList
//    {
//        public List<float[]> Magnitudes { get; set; }
//        public List<float[]> Phases { get; set; }
//    }

Samples can be reconstructed from Magnitude-Phase:

float[] reconstructed = stft.ReconstructMagnitudePhase(magPhase);

Window functions

13 well-known windowing functions are implemented in NWaves:

  • Rectangular
  • Triangular
  • Hamming
  • Blackman
  • Hann
  • Gaussian
  • Kaiser
  • Kaiser-Bessel derived
  • Bartlett-Hann
  • Lanczos
  • Power of sine
  • Flat-top
  • Liftering

Window samples can be generated one of two ways:

// 1)
var blackmanWindow = Window.Blackman(512);

// 2)
var hannWindow = Window.OfType(WindowTypes.Hann, 512);

WindowTypes is enum for window functions. Many algorithms in NWaves can accept parameters of this type.

All windows require length parameter, and some of them have additional parameters:

var kaiserWindow1 = Window.Kaiser(512);           // by default, alpha = 12
var kaiserWindow2 = Window.Kaiser(512, alpha: 9);
var kaiserWindow3 = Window.OfType(WindowTypes.Kaiser, 512, 15);

var kbdWindow1 = Window.Kbd(512);                 // by default, alpha = 4
var kbdWindow2 = Window.Kbd(512, alpha: 9);
var kbdWindow3 = Window.OfType(WindowTypes.Kbd, 512, 15);

var sineWindow1 = Window.PowerOfSine(512);        // by default, alpha = 1.5
var sineWindow2 = Window.PowerOfSine(512, alpha: 0.5);
var sineWindow3 = Window.OfType(WindowTypes.PowerOfSine, 512, 2);

var lifterWindow1 = Window.Liftering(512);        // by default, l = 22
var lifterWindow2 = Window.Liftering(512, l: 16);
var lifterWindow3 = Window.OfType(WindowTypes.Liftering, 512, 28);

There are several extension methods for applying windows:

var signal = new DiscreteSignal(8000, 128, 1f);
var samples = new float[128];

var window = Window.Blackman(128);

signal.ApplyWindow(window);    // 1) window samples are pre-computed
samples.ApplyWindow(samples);  // 2) window samples are pre-computed

signal.ApplyWindow(WindowTypes.Blackman);   // 3) window samples will be created each time
samples.ApplyWindow(WindowTypes.Blackman);  // 4) window samples will be created each time

Fast Wavelet Transform

Doing FWT is easy as this:

var input = new float[192];   // fill it with values
var output = new float[192];

var fwt = new Fwt(192, new Wavelet("db5"));

fwt.Direct(input, output);
fwt.Inverse(output, input);

FWT is compliant with pyWavelets/MATLAB given that signal extension mode is set there to 'periodization'/'ppd'.

FWT constructor accepts two parameters: FWT size and the Wavelet object. The size should not be necessarily a power of 2. However, currently, it must guarantee that at the given level of wavelet decomposition, the block of samples will have even size. Otherwise, no exception will be thrown but the results will be incorrect. The power-of-2 size guarantees that arrays will always have even size.

var fwt = new Fwt(108, new Wavelet("db5"));

var maxLevel = fwt.MaxLevel(108);  // maxLevel = 3
                                   // maxLevel is auto-computed and used for Fwt by default

// 108 samples ->   54 samples + 54 samples

fwt.Direct(input, output, 1);     // only first level, OK
fwt.Inverse(output, input, 1);    // only first level, OK

// ============================ or ===================================

// level1:  108 samples (even) ->   54 samples + 54 samples
// level2:  54 samples  (even) ->   27 samples + 27 samples

fwt.Direct(input, output, 2);     // only two levels, OK
fwt.Inverse(output, input, 2);    // only two levels, OK

// ============================ or ===================================

// level1:  108 samples (even) ->   54 samples + 54 samples
// level2:  54 samples  (even) ->   27 samples + 27 samples
// level3:  27 samples  (odd)  ->   13 samples + 14 samples

fwt.Direct(input, output, 3);   // at the third level results are incorrect
fwt.Inverse(output, input, 3);  // at the third level results are incorrect

Wavelets are represented with class Wavelet:

var w = new Wavelet(WaveletFamily.Symlet, 7);
// or
var w = new Wavelet("sym7");

var name = w.Name;     // "sym7"
var length = w.Length; // 14
var loD = w.LoD;       // { 0,002681815, -0,001047385, -0,0126363, ... }
var hiD = w.HiD;       // { ... }
var loR = w.LoR;       // { ... }
var hiR = w.HiR;       // { ... }

Supported wavelets are:

  • Haar ("haar")
  • Daubechies ("db1", "db2", "db3", ..., "db20")
  • Symlet ("sym2", "sym3", "sym4", ..., "sym20")
  • Coiflet ("coif1", "coif2", "coif3", "coif4", "coif5")

DCT transformers

There are several different well-known DCTs. In NWaves DCTs of type I-IV are implemented. They are fully compliant with scipy/MATLAB equivalents. DCT classes implement three methods of IDct interface:

  • Direct
  • Inverse
  • DirectNorm (normalized DCT)
var input = new float[20];  // fill it with values
var output = new float[12]; // take only 12 compressed values

var dct1 = new Dct1(20);    // type I
dct1.Direct(input, output);
dct1.DirectNorm(input, output);
dct1.Inverse(output, input);

var dct2 = new Dct2(20);   // type II
dct2.Direct(input, output);
dct2.DirectNorm(input, output);
dct2.Inverse(output, input);

var dct3 = new Dct3(20);   // type III
dct3.Direct(input, output);
dct3.DirectNorm(input, output);
dct3.Inverse(output, input);

var dct4 = new Dct4(20);   // type IV
dct4.Direct(input, output);
dct4.DirectNorm(input, output);
dct4.Inverse(output, input);

For large arrays DCT can be computed more efficiently via FFT.

Hartley transformer

var dht = new HartleyTransform(512);

dht.Direct(re);    // in-place
dht.Inverse(re);   // in-place

Mellin transformer

Mellin transformer has some specific features:

  • sizes of both input and output array must be specified (and (1) they may differ; (2) input size can be arbitrary, output size must be power of two)
  • currently, only Direct method is implemented

The following code demonstrates that mellin spectra of a signal and its time-stretched version are (almost) identical:

var stretch = 1.5;

int length = 256;
int stretchedLength = (int)(stretch * length);

// build sinusoid
var sin1 = new SineBuilder()
                   .SetParameter("freq", 5/*Hz*/)
                   .SampledAt(length)
                   .OfLength(length)
                   .Build();

// build the same sinusoid but time-stretched
var sin2 = new SineBuilder()
                   .SetParameter("freq", 5/*Hz*/)
                   .SampledAt(stretchedLength)
                   .OfLength(stretchedLength)
                   .Build();

var outputRe = new float[128];
var outputIm = new float[128];

// compute Mellin power spectrum of sinusoid
var mt1 = new MellinTransform(length, 128);
mt1.Direct(sin1.Samples, outputRe, outputIm);

var mellinSpectrum1 = outputRe.Take(64)
                              .Zip(outputIm, (re, im) => re * re + im * im)
                              .ToArray();

// compute Mellin power spectrum of time-stretched sinusoid
var mt2 = new MellinTransform(stretchedLength, 128);
mt2.Direct(sin2.Samples, outputRe, outputIm);

var mellinSpectrum2 = outputRe.Take(64)
                              .Zip(outputIm, (re, im) => re * re + im * im)
                              .ToArray();

// compute Mean Absolute Error between two spectra:
var err = mellinSpectrum1.Zip(mellinSpectrum2, (s1, s2) => Math.Abs(s1 - s2)).Sum();

Cepstral transformer:

var ct = new CepstralTransform(20, fftSize: 512, logBase: 10);

var delay = ct.Direct(signal, cceps);  // direct complex cepstrum:  signal -> cceps
ct.Inverse(cceps, signal, delay);      // inverse complex cepstrum: cceps -> signal

ct.RealCepstrum(signal, rceps);        // real cepstrum:  signal -> rceps
ct.PowerCepstrum(signal, pceps);       // power cepstrum: signal -> pceps
ct.PhaseCepstrum(signal, hceps);       // phase cepstrum: signal -> hceps

The first parameter in the constructor of CepstralTransform is the number of cepstral coefficients. The second (optional) parameter is the size of FFT that will be used to compute cepstrum. The third (optional) parameter is the base of logarithm (by default it's Math.E, i.e. natural logarithm is used for computations).

The results are compliant with output of MATLAB functions cceps/icceps and rceps.

Complex cepstrum is calculated according to formula:

img

Phase unwrapping is quite tricky and differs from standard MathUtils.Unwrap. Also, the delay parameter arises from this. Read more about it here.

Real cepstrum is computed similarly, with the term phase_unwrapped discarded.

Power cepstrum is calculated from real cepstrum as: power = 4 * real ^ 2.

Phase cepstrum is calculated from complex cepstrum as: phase = (complex_cepstrum - reversed_complex_cepstrum) ^ 2.

Hilbert transformer:

var ht = new HilbertTransform(1024);
var result = ht.Direct(doubleSamples);

HilbertTransform class also provides method for computing complex analytic signal. Thus, previous line is equivalent to:

var result = ht.AnalyticSignal(doubleSamples).Imag;

By default HilbertTransformer works with double precision. This code is for floats:

var ht = new HilbertTransform(1024, doublePrecision: false);
var result = ht.AnalyticSignal(floatSamples).Item2;

Goertzel algorithm

Example of computing only one term (#2) of 8-size FFT using Goertzel algorithm:

float[] array = { 1, 2, 3, 4, 5, 6, 7, 8 };

Complex secondFftTerm = new Goertzel(8).Direct(array, 2);

Efficient spectral processing

In previous cases the result of each transform was a newly created object of the DiscreteSignal class. If the sequence of blocks must be processed then it's better to work with reusable arrays in memory (all output and intermediate results will be stored in reusable arrays, i.e. there will be no unnecessary memory allocations):

var fft = new RealFft(1024);

var spectrum = new float[1024];
var cepstrum = new float[20];

fft.PowerSpectrum(signal[1000, 2024].Samples, spectrum);
// do something with spectrum

fft.PowerSpectrum(signal[2024, 3048].Samples, spectrum);
// do something with spectrum

fft.PowerSpectrum(signal[3048, 4072].Samples, spectrum);
// do something with spectrum

ct.Direct(signal[5000, 5512].Samples, cepstrum)
// do something with cepstrum

//...

Hence, the following code is the best for spectral analysis and widely used in many parts of NWaves library:

// for extension method FastCopyTo() 
using NWaves.Utils;

var fftSize = 512;

var fft = new RealFft(fftSize);
var window = Window.Blackman(fftSize);
var block = new float[fftSize];
var spectrum = new float[fftSize / 2 + 1];

// loop:

for (var pos = 0; pos + fftSize < signal.Length; pos += hopSize)
{
    // copy 512 samples to 'block' starting from index 'pos'
    signal.Samples.FastCopyTo(block, 512, pos);

    block.ApplyWindow(window);

    fft.MagnitudeSpectrum(block, spectrum);

    // post-process/visualize spectrum
}

64-bit FFT

In filter design and analysis we're dealing with doubles (64 bit). Therefore, 64-bit version of complex FFT was added to library. It can be used just like its 32-bit counterpart, however, only Direct() and Inverse methods are available:

// for extension method ToDoubles() 
using NWaves.Utils;

var fft = new Fft64(512);

double[] real = signal.First(512).Samples.ToDoubles();
double[] imag = new double [512];

// in-place FFT
fft.Direct(real, imag);

// ...do something with real and imaginary parts of the spectrum...

// in-place IFFT
fft.Inverse(real, imag);
Clone this wiki locally