{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":70203,"databundleVersionId":8068726,"sourceType":"competition"}],"dockerImageVersionId":30673,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"import numpy as np\nimport matplotlib.pyplot as plt\nimport matplotlib.cm as cm\nfrom scipy.fft import fft, ifft\nfrom IPython.display import Audio\n\ndef zoom_slice(t, slice_size):\n    center = len(t) // 2\n    return slice(center - slice_size, center + slice_size)","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2024-04-07T00:31:13.695027Z","iopub.execute_input":"2024-04-07T00:31:13.695847Z","iopub.status.idle":"2024-04-07T00:31:13.703352Z","shell.execute_reply.started":"2024-04-07T00:31:13.695808Z","shell.execute_reply":"2024-04-07T00:31:13.702310Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Low Pass Filter - Frequency Domain Multiplication Approach\nA low pass audio filter removes frequencies above some cut off point, and keeps the frequencies below that cut off. \nFirst of all we'll need an audio signal to work with. We can construct a simple signal using a sine wave. \n\nHere's a sine wave oscillating at 55 HZ (A1 on piano) over 1 second.","metadata":{}},{"cell_type":"code","source":"def plot_simple_sine_wave(frequency, duration, sample_rate):\n    t = np.arange(0, duration, 1/sample_rate)\n    sine_wave =  0.5 * np.sin(2 * np.pi * frequency * t)\n\n    plt.figure(figsize=(20, 5))\n    plt.plot(t, sine_wave)\n    plt.xlabel('Time (s)')\n    plt.ylabel('Amplitude')\n    plt.title(f\"Sine wave at {frequency} Hz\")\n    plt.show()\n\nplot_simple_sine_wave(55, 1, 220)\n","metadata":{"execution":{"iopub.status.busy":"2024-04-07T00:31:13.705681Z","iopub.execute_input":"2024-04-07T00:31:13.706375Z","iopub.status.idle":"2024-04-07T00:31:14.037383Z","shell.execute_reply.started":"2024-04-07T00:31:13.706315Z","shell.execute_reply":"2024-04-07T00:31:14.036440Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's take a listen to the sine wave","metadata":{}},{"cell_type":"code","source":"\ndef play_simple_sine_wave(freq, sample_rate, duration, amplitude=0.1):\n    t = np.linspace(0, duration, int(sample_rate * duration), endpoint=False)\n    y = amplitude * np.sin(2 * np.pi * freq * t)\n    return Audio(y, rate=sample_rate, normalize=False)\n\nplay_simple_sine_wave(55, 44100, 1, 0.1)","metadata":{"execution":{"iopub.status.busy":"2024-04-07T00:31:14.038896Z","iopub.execute_input":"2024-04-07T00:31:14.040019Z","iopub.status.idle":"2024-04-07T00:31:14.053595Z","shell.execute_reply.started":"2024-04-07T00:31:14.039977Z","shell.execute_reply":"2024-04-07T00:31:14.052281Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"To demonstrate a low pass filter we'll need a signal with multiple frequencies.\nHere's a sine wave with frequencies 55 HZ, 110 HZ, 220 HZ, 440 HZ, and 880 HZ","metadata":{}},{"cell_type":"code","source":"def multi_frequency_sine_wave(frequencies, duration, sample_rate, amplitude=0.1):\n    t = np.arange(0, duration, 1/sample_rate)\n    nf = len(frequencies)\n    sine_wave = np.zeros((nf, len(t)))\n    for i, f in enumerate(frequencies):\n        sine_wave[i] = (1 / nf) * np.sin(2 * np.pi * f * t)\n    return amplitude * sine_wave\n\ndef plot_multi_frequency_sine_wave(frequencies, duration, sample_rate):\n    sine_wave = multi_frequency_sine_wave(frequencies, duration, sample_rate)\n\n    fig, axs = plt.subplots(2, 1, figsize=(20, 10))\n    t = np.arange(0, duration, 1/sample_rate)\n    axs[0].set_title(f\"Signal of sine waves at {frequencies} Hz\")\n    axs[0].plot(t, sine_wave.T)\n\n    slice_size = 110\n    axs[1].set_title(f\"First {slice_size} samples\")\n    axs[1].plot(t[:slice_size], sine_wave.T[:slice_size])\n    axs[0].set_xlabel('Time (s)')\n    plt.show()\n    \n\ndef play_multi_frequency_sine_wave(frequencies, duration, sample_rate, amplitude=0.01):\n    sine_wave = multi_frequency_sine_wave(frequencies, duration, sample_rate, amplitude)\n    return Audio(sine_wave, rate=sample_rate, normalize=False)\n\nplot_multi_frequency_sine_wave([55, 110, 220, 440, 880], 1, 880 * 4)\nplay_multi_frequency_sine_wave([55, 110, 220, 440, 880], 2, 44100, 0.025)","metadata":{"execution":{"iopub.status.busy":"2024-04-07T00:31:14.057312Z","iopub.execute_input":"2024-04-07T00:31:14.057761Z","iopub.status.idle":"2024-04-07T00:31:15.011840Z","shell.execute_reply.started":"2024-04-07T00:31:14.057713Z","shell.execute_reply":"2024-04-07T00:31:15.010785Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Usually we will not know how a signal has been constructed so we'll have to use the discrete fourier transform to deconstruct the signal into it's constitute sine waves.\n\nThe fourier transform, transforms a signal from the time domain to the frequency domain using integration. The parts of this integral are sine and cosine waves. In the discrete case this integral becomes a sum. \n\n### Definition Of The Discrete Fourier Transform\n$$\nF(k) = \\sum_{n=0}^{N-1} f(n) \\cdot e^{-\\frac{2\\pi i}{N}kn}\n$$\nhttps://en.wikipedia.org/wiki/Discrete_Fourier_transform\n\nBoth the continuous and discrete fourier transforms are usually written in the complex exponential form. In other words it uses Euler's formula shown below to connect the trigonometric functions sine and cosine to the complex numbers. I recommend checking out Steven Bruton's excellent video on Euler's formula https://youtu.be/Rp-smPZLESc?si=Ru-hL4PaFcrZ5DeG\n\n### Definition Of Euler's Formula\n$$\ne^{ix} = \\cos(x) + i\\sin(x)\n$$\nhttps://en.wikipedia.org/wiki/Euler%27s_formula","metadata":{}},{"cell_type":"code","source":"class SignalGenerator:\n    def __init__(self, fs, duration, base_freq):\n        self.fs = fs\n        self.duration = duration\n        self.base_freq = base_freq\n        self.t = np.linspace(0, duration, int(fs*duration), endpoint=False)\n        \n    def generate_signal(self, num_harmonics):\n        self.signal = np.sum([np.sin(2*np.pi*(self.base_freq*2**i) * self.t) for i in range(num_harmonics)], axis=0)\n        return self.signal, self.t, self.fs\n    \n    def plot(self, zoom_slice):\n        signal_fft = fft(self.signal)\n        freq = np.fft.fftfreq(len(self.t), d=1/self.fs)\n\n        plt.figure(figsize=(16, 4))\n        plt.subplot(1, 2, 1)\n        plt.plot(self.t[zoom_slice], self.signal[zoom_slice])\n        plt.title('Original Signal')\n        plt.xlabel('Time [s]')\n        plt.ylabel('Amplitude')\n\n        plt.subplot(1, 2, 2)\n        plt.plot(freq, np.abs(signal_fft))\n        plt.title('Signal in Frequency Domain')\n        plt.xlabel('Frequency [Hz]')\n        plt.ylabel('Magnitude')\n\n        plt.tight_layout()\n        plt.show()\n        \n    \nsg = SignalGenerator(fs=8000, duration=1, base_freq=55)\nsignal, signal_t, signal_fs = sg.generate_signal(5)\nsg.plot(slice(0,100))\n","metadata":{"execution":{"iopub.status.busy":"2024-04-07T00:31:15.013090Z","iopub.execute_input":"2024-04-07T00:31:15.013388Z","iopub.status.idle":"2024-04-07T00:31:15.571252Z","shell.execute_reply.started":"2024-04-07T00:31:15.013363Z","shell.execute_reply":"2024-04-07T00:31:15.570215Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"So now we can see which frequencies the signal is made up of. To create a low pass filter we can remove some of the high frequencies. We can do this my multiplying the signals we want to remove by 0 and multiple the ones we want to keep by 1, in other words we multiple by a rectangle function. ","metadata":{}},{"cell_type":"code","source":"cutoff = 441\nsignal_fft = fft(signal)\nfreq = np.fft.fftfreq(len(signal_t), d=1/signal_fs)\nrect_filter = np.array([0 if abs(f) >= cutoff else 1 for f in freq])\nfiltered_fft = signal_fft * rect_filter\n\nplt.figure(figsize=(16, 4))\nplt.subplot(1, 2, 1)\nplt.plot(freq, rect_filter)\nplt.title('Rectangle Function In Frequency Domain')\nplt.xlabel('Frequency [Hz]')\nplt.ylabel('Magnitude')\n\nplt.subplot(1, 2, 2)\nplt.plot(freq, np.abs(filtered_fft))\nplt.title('Filtered Signal in Frequency Domain')\nplt.xlabel('Frequency [Hz]')\nplt.ylabel('Magnitude')\n\nplt.tight_layout()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-04-07T00:31:15.572637Z","iopub.execute_input":"2024-04-07T00:31:15.572949Z","iopub.status.idle":"2024-04-07T00:31:16.138575Z","shell.execute_reply.started":"2024-04-07T00:31:15.572924Z","shell.execute_reply":"2024-04-07T00:31:16.137431Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Ok so we have removed some of the high frequencies but now we need to get our signal back. We can get the signal back using the inverse fourier transform. ","metadata":{}},{"cell_type":"code","source":"# Inverse Fourier Transform\nfiltered_signal = ifft(filtered_fft)\n\nplt.figure(figsize=(16, 4))\n\nplt.subplot(1, 2, 1)\nplt.plot(signal_t[:100], filtered_signal.real[:100])  \nplt.title('Filtered Signal in Time Domain (Low Pass Filter Applied)')\nplt.xlabel('Time [s]')\nplt.ylabel('Amplitude')\n\nplt.subplot(1, 2, 2)\nplt.plot(signal_t[:100], signal[:100])\nplt.title('Original Signal')\nplt.xlabel('Time [s]')\nplt.ylabel('Amplitude')\n\nplt.tight_layout()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-04-07T00:31:16.139717Z","iopub.execute_input":"2024-04-07T00:31:16.140009Z","iopub.status.idle":"2024-04-07T00:31:16.659073Z","shell.execute_reply.started":"2024-04-07T00:31:16.139984Z","shell.execute_reply":"2024-04-07T00:31:16.658026Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"And there we have it. A low pass audio filter!\n\nHere's the playback of the filtered and original signal. Notice the fitlered signal does not contain the high frequency tones","metadata":{}},{"cell_type":"code","source":"def play_filtered_signal(signal, sample_rate):\n    normalized_signal = 0.01 * (signal / np.max(np.abs(signal)))\n    return Audio(normalized_signal, rate=sample_rate, normalize=False)\n\nplay_filtered_signal(filtered_signal, signal_fs)","metadata":{"execution":{"iopub.status.busy":"2024-04-07T00:31:16.660523Z","iopub.execute_input":"2024-04-07T00:31:16.660847Z","iopub.status.idle":"2024-04-07T00:31:16.669927Z","shell.execute_reply.started":"2024-04-07T00:31:16.660820Z","shell.execute_reply":"2024-04-07T00:31:16.668988Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def play_original_signal(signal, sample_rate):\n    return Audio(0.005 * signal, rate=sample_rate, normalize=False)\n\nplay_original_signal(signal, signal_fs)","metadata":{"execution":{"iopub.status.busy":"2024-04-07T00:31:16.671205Z","iopub.execute_input":"2024-04-07T00:31:16.671540Z","iopub.status.idle":"2024-04-07T00:31:16.681839Z","shell.execute_reply.started":"2024-04-07T00:31:16.671514Z","shell.execute_reply":"2024-04-07T00:31:16.680651Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Low Pass Filter - Time Domain Convolution Approach\nIn practice, especially for real time audio processing filtering is usually done via convolution in the time domain. The convolution theorem tells us that convolution of signals in the time domain corresponds to taking the Fourier transform of each signal and multipling them in the Fourier domain, then taking the inverse Fourier transform of the result.\n\nAbove we demonstrated filtering via multiplication in the Fourier domain. If we want to get the same result in the time domain using convolution we could take the inverse Fourier transform of our rectangle function and convolve the result with our signal in the time domain. But we can take a shortcut because the inverse fourier transform of the rect function happens to be the sinc(x) function!\n\n### Definition Of The Discrete Convolution\nNotice that like the Fourier transform, convolution is also an integral but it maps it back to the same domain, in this case time domain -> time domaine, whereas the Fourier transform maps from time domain -> frequency domain. \n$$\n(f * g)[n] = \\sum_{k=-\\infty}^{\\infty} f[k] \\cdot g[n - k]\n$$\nhttps://en.wikipedia.org/wiki/Convolution\n\n### The Convolution Theorem\nIf we convolve f with g and take the Fourier transform, we get the same result as taking the Fourier transform of f, taking the Fourier transform of g and multiplying the results.\nIn other words, convolution in the time domain corresponds to multiplication in the Fourier domain.\n$$\n\\mathcal{F}\\{f * g\\} = \\mathcal{F}\\{f\\} \\cdot \\mathcal{F}\\{g\\}\n$$\nhttps://en.wikipedia.org/wiki/Convolution_theorem\n\n### Rect and Sinc\nThe Fourier transform of the sinc function gives us the rectangle function.\n$$\n\\mathcal{F}\\{\\text{sinc}(t)\\} = \\text{rect}(\\omega)\n$$\n\nThe Fourier transform of the rectangle function gives us the sinc function.\n$$\n\\mathcal{F}\\{\\text{rect}(\\omega)\\} = \\text{sinc}(t)\n$$\nhttps://en.wikipedia.org/wiki/Sinc_function","metadata":{}},{"cell_type":"code","source":"cutoff = 220 \nN = 4096 \nfreq = np.fft.fftfreq(N)\n\n# Construct the rectangle function\nrect_filter = np.array([0 if abs(f * N) >= cutoff else 1 for f in freq])\n\n# Perform the inverse FFT\ntime_signal = np.fft.ifft(rect_filter)\n\n# Shift the time-domain signal to center it\ntime_signal_shifted = np.fft.ifftshift(time_signal)\n\n# Normalize the signal to have a peak amplitude of 1\ntime_signal_normalized = time_signal_shifted / np.max(np.abs(time_signal_shifted))\n\nt = np.arange(-N//2, N//2)\n\n# Scale the sinc function by the ratio of the cutoff frequency to half the number of points in the FFT\nvanilla_sinc_function = np.sinc(cutoff * t / (N // 2))\n\nplt.figure(figsize=(16, 4))\nplt.subplot(1, 3, 1)\nplt.plot(freq, rect_filter)\nplt.title('Rect Function in Frequency Domain')\nplt.xlabel('Frequency')\nplt.ylabel('Amplitude')\n\nplt.subplot(1, 3, 2)\nplt.plot(t[zoom_slice(t, 100)], np.real(time_signal_normalized)[zoom_slice(t, 100)])\nplt.title('Sinc Function From Inverse Fourier Transform of Rect Function')\nplt.xlabel('Time')\nplt.ylabel('Amplitude')\n\nplt.subplot(1, 3, 3)\nplt.plot(t[zoom_slice(t, 100)], vanilla_sinc_function[zoom_slice(t, 100)])\nplt.title('Vanilla Sinc Function')\nplt.xlabel('Time')\nplt.ylabel('Amplitude')\nplt.show()\n\n\nplt.tight_layout()\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2024-04-07T00:31:16.684893Z","iopub.execute_input":"2024-04-07T00:31:16.685221Z","iopub.status.idle":"2024-04-07T00:31:17.422773Z","shell.execute_reply.started":"2024-04-07T00:31:16.685193Z","shell.execute_reply":"2024-04-07T00:31:17.421578Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class ConvolutionLowPassFilter:\n    def __init__(self, cutoff, N):\n        self.cutoff = cutoff\n        self.N = N\n        \n    def apply_filter(self, signal):\n        t_sinc = np.arange(-N//2, N//2)\n\n        # Scale the sinc function by the ratio of the cutoff frequency to half the number of points in the FFT\n        sinc_function = np.sinc(cutoff * t_sinc / (N // 2))\n\n        # Convolve the signal with the sinc function\n        filtered_signal = np.convolve(signal, sinc_function, mode='same')\n\n        return filtered_signal\n\nlow_pass_filter = ConvolutionLowPassFilter(cutoff=220, N=4096)\n\n# Convolve the signal with the sinc function\nconvolution_filtered_signal = low_pass_filter.apply_filter(signal)\n\nplt.figure(figsize=(16, 4))\nplt.subplot(1, 2, 1)\nplt.plot(signal_t[zoom_slice(t, 110)], convolution_filtered_signal[zoom_slice(t, 110)], label='Filtered Signal')\nplt.title('Low Pass Filter Applied via Convolving Sinc Function with Signal')\nplt.xlabel('Time')\nplt.ylabel('Amplitude')\nplt.legend()\n\nplt.subplot(1, 2, 2)\nplt.plot(signal_t[zoom_slice(t, 110)], signal[zoom_slice(t, 110)], label='Original Signal')\nplt.title('Original Signal')\nplt.xlabel('Time')\nplt.ylabel('Amplitude')\nplt.legend()\n\nplt.tight_layout()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-04-07T00:31:17.424094Z","iopub.execute_input":"2024-04-07T00:31:17.424624Z","iopub.status.idle":"2024-04-07T00:31:17.984448Z","shell.execute_reply.started":"2024-04-07T00:31:17.424594Z","shell.execute_reply":"2024-04-07T00:31:17.983569Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"play_filtered_signal(convolution_filtered_signal, signal_fs)","metadata":{"execution":{"iopub.status.busy":"2024-04-07T00:31:17.985539Z","iopub.execute_input":"2024-04-07T00:31:17.986145Z","iopub.status.idle":"2024-04-07T00:31:17.991858Z","shell.execute_reply.started":"2024-04-07T00:31:17.986117Z","shell.execute_reply":"2024-04-07T00:31:17.991146Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"play_original_signal(signal, signal_fs)","metadata":{"execution":{"iopub.status.busy":"2024-04-07T00:31:17.992795Z","iopub.execute_input":"2024-04-07T00:31:17.993466Z","iopub.status.idle":"2024-04-07T00:31:18.005302Z","shell.execute_reply.started":"2024-04-07T00:31:17.993437Z","shell.execute_reply":"2024-04-07T00:31:18.004144Z"},"trusted":true},"execution_count":null,"outputs":[]}]}