drummyfish/less_retarded_wiki
5
1
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity

Update

Browse source
This commit is contained in:
Miloslav Ciz 2023-12-02 16:47:02 +01:00
parent 3ca4661a3d
commit d5b4366475
10 changed files with 199 additions and 10 deletions
Download patch file Download diff file Expand all files Collapse all files

186
fourier_transform.md
View file

@ -8,7 +8,7 @@ For newcomers FT is typically not easy to understand, it takes time to wrap one'
FT is actually just one of many so called **[integral transforms](integral_transform.md)** that are all quite similar -- they always transform the signal to some other domain and back, they use similar equation but usually use a different kind of [function](function.md). Other integral transforms are for example **[discrete cosine transformation](dct.md)** (DCT) or **[wavelet transform](wavelet_transform.md)**. DCT is actually a bit simpler than FT, so if you are having hard time with FT, go check out DCT.
**If you know [linear algebra](linear_algebra.md), this may help you understand what (D)FT really does:** Imagine the signal we work with is a POINT (we can also say a [vector](vector.md)) in many [dimensional](dimension.md) space; if for example we have a recorded sound that has 1000 samples, it is really a 1000 dimensional vector, a point in 1000 dimensional space, expressed as an "array" of 1000 numbers (vector components). A short note: since we consider a finite number of discrete samples here, we are actually dealing with what's called DISCRETE FT here, not the "vanilla" FT, but for now let's not diverge. (D)FT does nothing more than transforming from one vector [basis](basis.md) ("system of coordinates", "frame of reference") to another basis; i.e. by default the signal is expressed in time domain (our usual vector basis), the numbers in the sound "array" are such because we are viewing them from the time "frame of reference" -- (D)FT will NOT do anything with to the signal itself (it is a vector/point in space, which will stay where it is, the recorded sound itself will not change), it will merely express this same point/vector from a different "point of view"/"frame of reference" (set of basis vectors) -- that of frequencies. That's basically how all the integral transforms work, they just have to ensure the basis they are transforming to is orthogonal (i.e. kind of legit, "usable") of course.
**If you know [linear algebra](linear_algebra.md), this may help you understand what (D)FT really does:** Imagine the signal we work with is a POINT (we can also say a [vector](vector.md)) in many [dimensional](dimension.md) space; if for example we have a recorded sound that has 1000 samples, it is really a 1000 dimensional vector, a point in 1000 dimensional space, expressed as an "array" of 1000 numbers (vector components). A short note: since we consider a finite number of discrete samples here, we are actually dealing with what's called DISCRETE FT here, not the "vanilla" FT, but for now let's not diverge. (D)FT does nothing more than transforming from one vector [basis](basis.md) ("system of coordinates", "frame of reference") to another basis; i.e. by default the signal is expressed in time domain (our usual vector basis), the numbers in the sound "array" are such because we are viewing them from the time "frame of reference" -- (D)FT will NOT do anything with to the signal itself (it is a vector/point in space, which will stay where it is, the recorded sound itself will not change), it will merely express this same point/vector from a different "point of view"/"frame of reference" (set of basis vectors) -- that of frequencies. That's basically how all the integral transforms work, they just have to ensure the basis they are transforming to is orthogonal (i.e. kind of legit, "usable") of course. In addition the FT equation is nothing complex, it literally just uses a **[dot product](dot_product.md)** of the whole input signal with each checked frequency wave to find out how similar the signal is to that particular frequency, as dot product simply says "how similar two vectors are" -- really, think about the equation and you will see it's really doing just that.
## Details
@ -22,4 +22,186 @@ First let's make clearer the whole terminology around FT:
- **Discrete Fourier Transform (DFT)** (not to be confused with DTFT!): Uses DFS to transform a FINITE DISCRETE signal to a FINITE DISCRETE spectrum (with the same period as the input) by simply "pretending" the finite input signal is actually repeating over and over and then, after the transform, only leaving in the first period of the result (since the rest is just repeating). **This is actually what programmers usually mean by Fourier Transform** because in computers we practically always only deal with finite discrete signals (i.e. [arrays](array.md) of data).
- **Fast Fourier Transform (FFT)**: Computes DFT (NOT FT!) that's faster than the [naive](naive.md) implementation, i.e. computing the equation that defines DFT as it's written has time complexity O(n^2) while FFT improves this to O(n * log(n)).
TODO: equations, more, code, pictures etc.
From now on we will implicitly be talking about DFT of a real function (we'll ignore the possibility of complex input), the most notable transform here.
The input to DFT is a real function, i.e. the time domain representation of the signal. The output is a complex valued function of frequency, i.e. the spectrum -- for each frequency it says a complex number whose magnitude and phase say the magnitude and phase of that frequency (a sine wave) in the signal (many programs will visualize just the magnitude part as that's usually the important thing, however keep in mind there is always also the phase part as well).
The general equations defining DFT and IDFT, for signal with *N* samples, are following
```
___ N - 1
\
DFT[k] = /__ x[n] * e^(-2 * i * pi * k * n / N)
n = 0
___ N - 1
\
IDFT[k] = 1/N * /__ x[n] * e^(2 * i * pi * k * n / N)
n = 0
```
OK, this is usually where every noob ragequits if he hasn't already because of all the [pi](pi.md)s and [e](e.md)s and just generally ununderstable mess of weird symbols etc. What the heck does this all mean? As said above, it's doing nothing else than [dot product](dot_product.md) or vectors really: one vector is the input signal and the other vectors are the individual frequencies (sine waves) we are trying to discover in the signal -- this looks so complicated because here we are actually viewing the general version for a possible [complex](complex_number.md) input signal, the *e to something* part is actually the above mentioned complex exponential, it is the exponential way of writing a complex number (see e.g. [Euler's identity](eulers_indentity.md)). Anyway, considering only real input signal, we can simplify this to a more programmer friendly form:
```
DFT:
init DFT_real and DFT_imag to 0s
for k = 0 to N - 1
for n = 0 to N - 1
angle = -2 * i * pi * k * n / N
DFT_real[k] += x[n] * cos(angle)
DFT_imag[k] += x[n] * sin(angle)
IDFT:
init data to 0s
for k = 0 to N - 1
for n = 0 to N - 1
angle = 2 * i * pi * k * n / N
data[k] += DFT_real[n] * cos(angle) - DFT_imag[n] * sin(angle)
data[k] /= N
```
**Example**: take a look at the following array of 8 kind of arbitrary values and what their DFT looks like:
```
# #
# # #
# # # #
# # # #
# # # # #
# # # # # #
data: 5.00 4.71 6.00 6.54 1.00 2.29 0.00 -0.54
DFT:
#
#
#
# # #
# # #
# # # # #
magn.: 25.00 12.74 1.00 7.33 1.00 7.33 1.00 12.74
phase: 0.00 -1.52 -1.57 -0.10 -3.14 0.10 1.57 1.52
-----
real: 25.00 0.70 0.00 7.30 -1.00 7.30 -0.00 0.70
imag.: 0.00 -12.72 -1.00 -0.72 -0.00 0.72 1.00 12.72
restored:
data: 5.00 4.71 6.00 6.54 1.00 2.29 0.00 -0.54
```
At the top we have the input data: notice the data kind of looks similar to a low-frequency sine wave, so the frequencies in the spectrum below are mostly low, but there's also some high frequency noise that's deforming the wave. For convenience here we show the spectrum values in both formats (magnitude/phase and real/imaginary part), but keep in mind it's just different formats of the same complex number values; for analysis we are mostly interested in the magnitude of the complex numbers as that shows as the amplitude of the frequency, i.e. the "amount" of the frequency in the signal. Here we notice the greatest peak is at frequency 0 -- this is a "constant" component, the lowest possible frequency that just represents a constant vertical offset of the signal (a constant number added to all samples); this component here is so big because our input signal doesn't really oscillate around the value 0 as it doesn't even go to negative values -- DFT sees this as our signal being shifted quite a lot "up". Frequencies 1 and 7 are the second biggest here: DFT is telling us the signal looks mostly like an addition of a sine wave with very low frequency and very high frequency (which it does), it doesn't see many middle value frequencies here. At the end we also see the original values computed back using IDFT, just to check everything is working as expected.
Here is the [C](c.md) code that generates the above, you may use it as a snippet and/or to play around with different inputs to see what their spectra look like (for "readability" we commit the sin of using [floating point](float.md) numbers here, implementation of DFT [without floats](fixed_point.md) is left as an exercise :]):
```
#include <stdio.h>
#include <math.h>
#define PI 3.141592
#define N (sizeof(data) / sizeof(double)) // size of input data
#define NUM_FORMAT "%6.2lf"
#define STR_FORMAT "%-10s"
#define DRAW_HEIGHT 6
double data[] = // enter input data here
{5.00, 4.71, 6.00, 6.54, 1.00, 2.29, 0.00, -0.54};
// output DFT:
double dftR[N]; // real part of DFT
double dftI[N]; // imaginary part of DFT
// just for printing:
double dftM[N]; // magnitude of DFT
double dftA[N]; // argument (angle/phase) of DFT
void printArray(double *array)
{
for (int i = 0; i < N; ++i)
printf(" " NUM_FORMAT,array[i]);
putchar('\n');
}
void drawArray(double *array, double scale)
{
for (int y = 0; y < DRAW_HEIGHT; ++y)
{
printf(" ");
for (int x = 0; x < N; ++x)
{
printf(" ");
putchar(((int) array[x] * scale) >= (DRAW_HEIGHT - y) ? '#' : ' ');
}
putchar('\n');
}
}
void printDft(void)
{
printf(STR_FORMAT," magn.:"); printArray(dftM);
printf(STR_FORMAT," phase:"); printArray(dftA);
puts(" -----");
printf(STR_FORMAT," real:"); printArray(dftR);
printf(STR_FORMAT," imag.:"); printArray(dftI);
}
void dft(void)
{
for (int i = 0; i < N; ++i)
{
dftR[i] = 0;
dftI[i] = 0;
}
for (int k = 0; k < N; ++k)
{
for (int n = 0; n < N; ++n)
{
double angle = (-2 * PI * k * n) / N;
dftR[k] += data[n] * cos(angle);
dftI[k] += data[n] * sin(angle);
}
// just for printing also precompute magnitudes and phases
dftM[k] = sqrt(dftR[k] * dftR[k] + dftI[k] * dftI[k]);
dftA[k] = atan2(dftI[k],dftR[k]);
}
}
void idft(void)
{
for (int i = 0; i < N; ++i)
data[i] = 0;
for (int k = 0; k < N; ++k)
{
for (int n = 0; n < N; ++n)
{
double angle = (2 * PI * k * n) / N;
data[k] += dftR[n] * cos(angle) - dftI[n] * sin(angle);
}
data[k] /= N;
}
}
int main(void)
{
drawArray(data,1);
printf(STR_FORMAT,"data:"); printArray(data);
puts("\nDFT:");
dft();
drawArray(dftM,0.25);
printDft();
idft();
puts("\nrestored:");
printf(STR_FORMAT,"data:"); printArray(data);
return 0;
}
```
TODO: pictures, 2D version