This gist is the implementation of Frequency-Division Multiplexing (FDM) to three audios with the same length and frequency sample

fdm.py

import wave
from playsound import playsound

from matplotlib import pyplot as plt
import numpy as np


def save_audio(x, fs, name):
    with wave.open(name, mode='wb') as wav:
        wav.setnchannels(1)  # Mono
        wav.setsampwidth(2)  # 16-bit
        wav.setframerate(fs)  # Sampling frequency

        wav.writeframes(x.astype('int16').tobytes())  # Write frames


def read_audio(audio):
    with wave.open(audio, mode='rb') as wav:
        frames = wav.readframes(wav.getnframes())  # Read all frames
        x = np.frombuffer(frames, dtype='int16')  # Convert from bytes to int16

        if wav.getnchannels() == 2:  # If stereo
            x = x[0::2]

        # x = x/max(x) * 1000  # Normalize and scale to 10000

        fs = wav.getframerate()  # Get sampling frequency
        dt = 1/fs  # Get sampling period

        t = np.arange(0, len(x)*dt, dt)  # Create time vector

    return x, t, dt


def shift(x, f, shift):
    x_shift = np.zeros_like(x)

    idx = np.where(f == shift)[0][0]
    x_shift = np.roll(x, idx)

    return x_shift


def bandpass(x, f, f_min, f_max):
    x_bandpass = np.zeros_like(x, dtype='complex')

    idx_min = np.where(f == f_min)[0][0]
    idx_max = np.where(f == f_max)[0][0]

    if f_min * f_max < 0:
        x_bandpass[idx_min:] = x[idx_min:]
        x_bandpass[:idx_max] = x[:idx_max]
    else:
        x_bandpass[idx_min:idx_max] = x[idx_min:idx_max]

    return x_bandpass


if __name__ == "__main__":
    audios = ["sine-440.wav", "square-440.wav", "triangle-440.wav"]
    shifts = [-15000, 0, 15000]

    X_m = None

    for i, audio in enumerate(audios):
        x, t, dt = read_audio(audio)
        # playsound(audio, block=True)

        X = np.fft.fft(x)  # Compute FFT
        f = np.fft.fftfreq(len(x), dt)  # Compute frequency vector for FFT

        X_shifted = shift(X, f, shift=shifts[i])  # Shift FFT

        if X_m is None:
            X_m = X_shifted.copy()
        else:
            X_m += X_shifted

        plt.figure(num=1)  # Plot audio time domain
        plt.subplot(len(audios), 1, i+1)
        plt.title(audio.replace("-440.wav", ""))
        plt.xlabel("Tiempo (s)")
        plt.ylabel("Amplitud")
        plt.plot(t, x)

        plt.figure(num=2)  # Plot audio frequency domain
        plt.subplot(len(audios), 1, i+1)
        plt.title(audio.replace("-440.wav", ""))
        plt.xlabel("Frecuencia (Hz)")
        plt.ylabel("Amplitud")
        plt.plot(f, np.abs(X))

    x_m = np.fft.ifft(X_m)  # Compute IFFT

    plt.figure(num=3)  # Plot audio time domain
    plt.subplot(2, 1, 1)
    plt.title("x_m")
    plt.xlabel("Tiempo (s)")
    plt.ylabel("Amplitud")
    plt.plot(t, x)

    plt.subplot(2, 1, 2)
    plt.title("x_m")
    plt.xlabel("Frecuencia (Hz)")
    plt.ylabel("Amplitud")
    plt.plot(f, np.abs(X_m))

    # Reconstruct audio
    for i, f_shift in enumerate(shifts):
        bandwidth = 10000
        X_filtered = bandpass(X_m, f, f_shift-bandwidth/2, f_shift+bandwidth/2)

        X_unshifted = shift(X_filtered, f, -f_shift)
        x_unshifted = np.fft.ifft(X_unshifted)

        name = audios[i].replace("-440.wav", "_reconstructed")
        save_audio(x_unshifted.real, 44100, name+".wav")

        plt.figure(num=4)  # Plot audio time domain
        plt.subplot(len(audios), 1, i+1)
        plt.title(name)
        plt.xlabel("Tiempo (s)")
        plt.ylabel("Amplitud")
        plt.plot(t, x_unshifted.real)

        plt.figure(num=5)  # Plot audio frequency domain
        plt.subplot(len(audios), 1, i+1)
        plt.title(name)
        plt.xlabel("Frecuencia (Hz)")
        plt.ylabel("Amplitud")
        plt.plot(f, np.abs(X_unshifted))

    plt.show()