Everything should be made as simple as possible, but not simpler. (Albert Einstein)

Thursday, December 26, 2013

Fast Fourier Transform in ChucK to Get Spectrum

We can convert time domain to frequency domain by using Fast Fourier Transform, that's to calculate the discrete fourier transformation. http://chuck.cs.princeton.edu/doc/language/uana.html

This is the basic of finding the spectrum of an input signal, in this case I use a sinusoidal signal as example.
https://github.com/flyingdisc/music-chuck/tree/master/FFTSpectrum
This will display frequencies and power (in polar form of the complex number) at the FFT bin indices.

// get FFT spectrum
// calculate frequency (approx)
// first bin is DC (freq = 0Hz)
// make sure to filter out the DC (signal - mean(signal) => signal)
// also make sure to applying band-pass filter at Nyquist frequency
// Nyquist frequency == sampling rate / 2

// Input signal
SinOsc g => FFT fft => blackhole;
// set samplingRate
second / samp => float samplingRate;
<<< "Sampling rate =", samplingRate >>>;

// FFT bin size
16 => fft.size;

// spectrum, first half bins, 0..N/2-1
// the rest is useless, it's only conjugate of the first half
complex spectral[fft.size()/2];

// a sample frequency of the sinusoidal input
5500 => g.freq;
    
// let fft.size samples pass
fft.size()::samp => now;
    
// take fast fourier transform
fft.upchuck();
// get the spectrum
fft.spectrum( spectral );

// display spectrum at Nth bin
<<< "Spectrum:" >>>;
for ( 0 => int i; i < fft.size()/2; i++ ) 
{
    i * samplingRate / fft.size() => float freq;
    <<< "bin:", i, "freq:", freq, "power:", spectral[i]$polar >>>;
}

The main frequency of signal has the highest magnitude of power, so we can finding it this way,

// get FFT spectrum
// calculate frequency (approx)
// first bin is DC (freq = 0Hz)
// make sure to filter out the DC (signal - mean(signal) => signal)
// also make sure to applying band-pass filter at Nyquist frequency
// Nyquist frequency == sampling rate / 2

// input signal
SinOsc g => FFT fft => blackhole;
// set samplingRate
second / samp => float samplingRate;
<<< "Sampling rate =", samplingRate >>>;

// FFT bin size
1024 => fft.size;

// spectrum, first half bins, 0..N/2-1
// the rest is useless, it's only conjugate of the first half
complex spectral[fft.size()/2];

// a sample frequency of the sinusoidal input
5500 => g.freq;
    
// let fft.size samples pass
fft.size()::samp => now;
    
// take fast fourier transform
fft.upchuck();
// get the spectrum
fft.spectrum( spectral );

// calculate highest power peak, 
// which is the main frequency
0 => float highestPower;
-1 => int highestIndex;
for ( 0 => int i; i < fft.size()/2; i++ ) 
{
    if ((spectral[i]$polar).mag > highestPower) 
    {
        (spectral[i]$polar).mag => highestPower;
        i => highestIndex;
    }
}
// calculate freq approximation
highestIndex * samplingRate / fft.size() => float freq0;
<<< "Frequency (approx):" >>>;
<<< "bin:", highestIndex, "freq:", freq0 >>>;
<<< "magnitude:", (spectral[highestIndex]$polar).mag, 
    "phase:", (spectral[highestIndex]$polar).phase >>>;

In the example, we use sinusoidal with frequency = 5500 Hz, the output of both codes above,

Sampling rate = 44100.000000 
"Spectrum:" : (string)
bin: 0 freq: 0.000000 power: %(0.0013,0.2437*pi) 
bin: 1 freq: 2756.250000 power: %(0.0018,-0.3305*pi) 
bin: 2 freq: 5512.500000 power: %(0.5004,0.5040*pi) 
bin: 3 freq: 8268.750000 power: %(0.0027,0.3952*pi) 
bin: 4 freq: 11025.000000 power: %(0.0015,0.3060*pi) 
bin: 5 freq: 13781.250000 power: %(0.0012,0.2251*pi) 
bin: 6 freq: 16537.500000 power: %(0.0010,0.1484*pi) 
bin: 7 freq: 19293.750000 power: %(0.0009,0.0737*pi) 

Sampling rate = 44100.000000 
"Frequency (approx):" : (string)
bin: 128 freq: 5512.500000 
magnitude: 0.433066 phase: 2.481109 

We see that the result is FFT bin with frequency 5512 Hz, for our sinusoidal 5500 Hz.

Note: 0th bin has DC component, that's when frequency = 0 Hz. And the power magnitude is RMS (root mean squared, although I'm not very sure about this).

I've not tried it yet for adc, but I think it would be the same, that's to replace SinOsc with adc, and removing the line with frequency setting.
Changed,
adc => FFT fft => blackhole;
Removed,
5500 => g.freq;

When applying adc as input, make sure that we remove (too high) DC component.
Also make sure to filter out all frequency components above Nyquist frequency (= sampling rate / 2).
I think this can be applied by using BPF (band pass filter) when needed. http://chuck.cs.princeton.edu/doc/program/ugen_full.html#BPF

Example, for audible sound [20Hz .. 20kHz], the sampling rate according to Nyquist criterion should be at least (2 x 20kHz) = 40kHz. Hence we must filter the input in the bandwidth of [20..20k] Hz, by using BPF (either by hardware or software).
In the example codes above, the sampling rate is 44.1kHz

Also for adc, we need also to change the FFT size as desired, e.g. to 1024 or 2048, the higher the size we get more smooth samples in frequency domain.
Also the duration,
fft.size()::samp => now;
need to be changed to cover the input signal duration.