{
"cells": [
{
"cell_type": "markdown",
"id": "4a93eae1-5c0e-4cf1-846a-621d0fcdd29e",
"metadata": {},
"source": [
"# Fast Fourier Transforms (FFT)\n",
"\n",
"Let $x_0, \\ldots, x_{n-1}$ be complex numbers. The [Discrete Fourier Transform (DFT)](https://en.wikipedia.org/wiki/Discrete-time_Fourier_transform) is defined by the formula:\n",
"\n",
"$$\n",
" X_k = \\sum_{m=0}^{n-1} x_m e^{-i2\\pi k m/n} \\qquad k = 0,\\ldots,n-1,\n",
"$$\n",
"\n",
"Let's define a stream, consisting of exactly 20 sine waves in 1024 samples:"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "554cb334-283a-40c1-89aa-5fe1fd41161c",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
]
},
"metadata": {
"": ""
},
"output_type": "display_data"
}
],
"source": [
"sines = Stream.SinOsc(samplerate*20/1024):sub(1, 1024)\n",
"sines:gnuplot()"
]
},
{
"cell_type": "markdown",
"id": "580e13a3-f910-472c-b887-76929aab1d87",
"metadata": {},
"source": [
"The bins of the corresponding frequency spectrum correspond to frequencies according to the following formula:\n",
"\n",
"$$\n",
" f_i = {i {f_{srate} \\over N}}\n",
"$$\n",
"\n",
"Plotting the frequency spectrum's magnitudes of `sines` should therefore give us a clear signal around bin 20 (Lua index 21).\n",
"\n",
"**NOTE:** This does not need windowing since the waves fit perfectly into 1024 sample buffer."
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "02aec3d0-4ae9-4bcb-a034-764418b7827f",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"{0, 0, 0, 1.173742634399e-14, 0, 0, 0, 3.9584680075211e-15, 0, 0, 0, 512, 0, 0, 0, 2.2315199565268e-14, 0, 0, 0, 8.6716166150709e-15, 0}\n"
]
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
]
},
"metadata": {
"": ""
},
"output_type": "display_data"
}
],
"source": [
"spectrum = tostream(magnitude(FFT(sines)))\n",
"print(spectrum:sub(10, 30))\n",
"spectrum:gnuplot()"
]
},
{
"cell_type": "markdown",
"id": "7ca7e0b3-1854-41b4-813f-016eb148387e",
"metadata": {},
"source": [
"Testing the inverse FFT - looks just like the original wave:"
]
},
{
"cell_type": "code",
"execution_count": 13,
"id": "e7fc0ea8-d510-41aa-815e-f29055f5b8e1",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
]
},
"metadata": {
"": ""
},
"output_type": "display_data"
}
],
"source": [
"IFFT(FFT(sines)):gnuplot()"
]
},
{
"cell_type": "markdown",
"id": "6c306c6e-dddc-4d8c-bac8-04f009c406df",
"metadata": {},
"source": [
"With 20.5 sines in the fragment, the plot is much better when applying a windowing function:"
]
},
{
"cell_type": "code",
"execution_count": 14,
"id": "03d8fa36-aad6-4de9-9866-1551c5010d25",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
]
},
"metadata": {
"": ""
},
"output_type": "display_data"
}
],
"source": [
"tostream(magnitude(FFT(Hamming(Stream.SinOsc(samplerate*20.5/1024):sub(1, 1024))))):gnuplot()"
]
},
{
"cell_type": "markdown",
"id": "c1e86db7-d3a6-4626-a71c-66bea1c04ae2",
"metadata": {},
"source": [
"Testing phase plotting. This is 20 sines but with phase 0.7 ($0.7 \\times 2\\pi$).\n",
"\n",
"**FIXME:** How exactly to interpret these numbers?"
]
},
{
"cell_type": "code",
"execution_count": 18,
"id": "0d8e9009-255c-4840-b930-408eb3c419ed",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.9390625, 0, 0, 0, 0, 0, 0, 0, 0, 0}\n"
]
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
]
},
"metadata": {
"": ""
},
"output_type": "display_data"
}
],
"source": [
"spectrum = tostream(phase(FFT(Stream.SinOsc(samplerate*20/1024, 0.7):sub(1, 1024))))\n",
"print(spectrum:sub(10, 30))\n",
"spectrum:gnuplot()"
]
},
{
"cell_type": "markdown",
"id": "12ed863d-99a4-4372-bb11-e0ce4b5013d3",
"metadata": {},
"source": [
"Let's try some naive real-time pitch shifting.\n",
"This however is [not easy to get right](https://www.reddit.com/r/DSP/comments/k6t24c/pitch_shifting_algorithm_in_frequency_domain/)!"
]
},
{
"cell_type": "code",
"execution_count": 20,
"id": "bd0e1d6d-ad05-46dc-9701-92846bc8a683",
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"WARNING: Buffer underrun detected\n"
]
}
],
"source": [
"-- FIXME: Use a different sound file!\n",
"robin = SndfileStream(\"tracks/robin-mono.wav\"):sub(sec(10.284), sec(17.466)):eval()\n",
"robin:play()\n",
"robin:FFT(1024):map(function(spectrum)\n",
" assert(#spectrum == 513)\n",
" -- NOTE: We cannot use Stream:resample() as it won't work with complex samples.\n",
" for i = 1, 512/2 do spectrum[i] = spectrum[i*2] end\n",
" for i = 512/2+1, 512 do spectrum[i] = 0i end\n",
" return spectrum\n",
"end):IFFT(1024):play()"
]
},
{
"cell_type": "markdown",
"id": "cdf40dfb-b7ac-42ae-a6c3-178c7c8db1ff",
"metadata": {},
"source": [
"Let's generate and plot the spectrum of a sine mixed with some noise:"
]
},
{
"cell_type": "code",
"execution_count": 21,
"id": "58845a87-ed29-45e8-85c7-ac4b24013e5c",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
]
},
"metadata": {
"": ""
},
"output_type": "display_data"
}
],
"source": [
"noisy = Stream.SinOsc(440):mix(NoiseStream, 0.4)\n",
"tostream(magnitude(FFT(Hamming(noisy:sub(1, 1024))))):mul(0.05):gnuplot()"
]
},
{
"cell_type": "markdown",
"id": "2141e582-948c-454e-834b-8148465e3869",
"metadata": {},
"source": [
"Try some very naive noise cancelling.\n",
"\n",
"**WARNING:** This plays indefinitely. Interrupt via Kernel → \"Interrupt Kernel\"."
]
},
{
"cell_type": "code",
"execution_count": 22,
"id": "b43ca02e-16b8-4738-a157-db6dc1bb5b22",
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"WARNING: Buffer underrun detected\n"
]
},
{
"ename": "n/a",
"evalue": "fft.lua:124: interrupted!",
"execution_count": 22,
"output_type": "error",
"traceback": [
"fft.lua:124: interrupted!",
"stack traceback:",
"\tfft.lua:124: in function 'FFT_internal'",
"\tfft.lua:119: in function 'FFT_internal'",
"\tfft.lua:119: in function 'FFT_internal'",
"\tfft.lua:136: in function 'fnc'",
"\tapplause.lua:1871: in function 'tick'",
"\tapplause.lua:1870: in function 'tick'",
"\tapplause.lua:1870: in function 'stream_tick'",
"\tapplause.lua:1634: in function 'tick'",
"\tapplause.lua:726: in function ",
"\t[C]: in function 'xpcall'",
"\tapplause.lua:683: in function ",
"\t[C]: in function 'xpcall'",
"\t...ies/ilua/env/lib/python3.8/site-packages/ilua/interp.lua:65: in function 'handle_execute'",
"\t...ies/ilua/env/lib/python3.8/site-packages/ilua/interp.lua:176: in main chunk",
"stack traceback:",
"\t[C]: in function 'error'",
"\tapplause.lua:707: in function ",
"\t[C]: in function 'xpcall'",
"\t...ies/ilua/env/lib/python3.8/site-packages/ilua/interp.lua:65: in function 'handle_execute'",
"\t...ies/ilua/env/lib/python3.8/site-packages/ilua/interp.lua:176: in main chunk"
]
}
],
"source": [
"noisy:FFT(1024, Hamming):map(function(spectrum)\n",
" assert(#spectrum == 513)\n",
" for i = 1, #spectrum do\n",
" if (spectrum[i].re^2 + spectrum[i].im^2)^.5*0.05 < 0.8 then spectrum[i] = 0i end\n",
" end\n",
" return spectrum\n",
"end):IFFT(1024):play()"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "fc38550e-835c-47a8-a792-7164109dafee",
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Lua",
"language": "lua",
"name": "lua"
},
"language_info": {
"file_extension": ".lua",
"mimetype": "text/x-lua",
"name": "lua",
"version": "n/a"
}
},
"nbformat": 4,
"nbformat_minor": 5
}