Copy code examples from mc-audio-effects

requirements.txt

@@ -1,4 +1,6 @@
 jupyter
 matplotlib
 numpy
+scipy
+sympy
 tabulate

scripts/adaa.py

@@ -0,0 +1,133 @@
+import matplotlib.pyplot as plt
+import numpy as np
+import scipy.signal as signal
+
+
+class ADAA1:
+    def __init__(self, f, AD1, TOL=1.0e-5):
+        self.TOL = TOL
+        self.f = f
+        self.AD1 = AD1
+
+    def process(self, x):
+        y = np.copy(x)
+        x1 = 0.0
+        for n, _ in enumerate(x):
+            if np.abs(x[n] - x1) < self.TOL:  # fallback
+                y[n] = self.f((x[n] + x1) / 2)
+            else:
+                y[n] = (self.AD1(x[n]) - self.AD1(x1)) / (x[n] - x1)
+            x1 = x[n]
+        return y
+
+
+class ADAA2:
+    def __init__(self, f, AD1, AD2, TOL=1.0e-5):
+        self.TOL = TOL
+        self.f = f
+        self.AD1 = AD1
+        self.AD2 = AD2
+
+    def process(self, x):
+        y = np.copy(x)
+
+        def calcD(x0, x1):
+            if np.abs(x0 - x1) < self.TOL:
+                return self.AD1((x0 + x1) / 2.0)
+            return (self.AD2(x0) - self.AD2(x1)) / (x0 - x1)
+
+        def fallback(x0, x2):
+            x_bar = (x0 + x2) / 2.0
+            delta = x_bar - x0
+
+            if delta < self.TOL:
+                return self.f((x_bar + x0) / 2.0)
+            return (2.0 / delta) * (self.AD1(x_bar) + (self.AD2(x0) - self.AD2(x_bar)) / delta)
+
+        x1 = 0.0
+        x2 = 0.0
+        for n, _ in enumerate(x):
+            if np.abs(x[n] - x1) < self.TOL:  # fallback
+                y[n] = fallback(x[n], x2)
+            else:
+                y[n] = (2.0 / (x[n] - x2)) * (calcD(x[n], x1) - calcD(x1, x2))
+            x2 = x1
+            x1 = x[n]
+        return y
+
+
+def plot_fft(x, fs, sm=1.0/24.0):
+    fft = 20 * np.log10(np.abs(np.fft.rfft(x) + 1.0e-9))
+    freqs = np.fft.rfftfreq(len(x), 1.0 / fs)
+    return freqs, fft
+
+
+def process_nonlin(fc, FS, nonlin, gain=10):
+    N = 200000
+    sin = np.sin(2 * np.pi * fc / FS * np.arange(N))
+    y = nonlin(gain * sin)
+    freqs, fft = plot_fft(y, FS)
+    return freqs, fft
+
+
+def signum(x):
+    return int(0 < x) - int(x < 0)
+
+
+def hardClip(x):
+    return np.where(x > 1, 1, np.where(x < -1, -1, x))
+
+
+def hardClipAD1(x):
+    return x * x / 2.0 if np.abs(x) < 1 else x * signum(x) - 0.5
+
+
+def hardClipAD2(x):
+    return x * x * x / 6.0 if np.abs(x) < 1 else ((x * x / 2.0) + (1.0 / 6.0)) * signum(x) - (x/2)
+
+
+def tanhClip(x):
+    return np.tanh(x)
+
+
+def tanhClipAD1(x):
+    return x - np.log(np.tanh(x) + 1)
+
+
+def softClip(x):
+    return -4*x**3/27 + x
+
+
+def softClipAD1(x):
+    return -x**4/27 + x**2/2
+
+
+FC = 1244.5
+FS = 44100
+OS = 2
+
+hardClip_ADAA1 = ADAA1(hardClip, hardClipAD1, 1.0e-5)
+hardClip_ADAA2 = ADAA2(hardClip, hardClipAD1, hardClipAD2, 1.0e-5)
+
+freqs_analog, fft_analog = process_nonlin(FC, FS*100, hardClip)
+freqs_alias, fft_alias = process_nonlin(FC, FS, hardClip)
+freqs_os, fft_os = process_nonlin(FC, FS*OS, hardClip)
+freqs_ad1, fft_ad1 = process_nonlin(FC, FS*OS, hardClip_ADAA1.process)
+freqs_ad2, fft_ad2 = process_nonlin(FC, FS*OS, hardClip_ADAA2.process)
+peak_idxs = signal.find_peaks(fft_analog, 65)[0]
+
+# plt.plot(freqs_analog, fft_analog, '--', c='red', label='Analog')
+# plt.plot(freqs_alias, fft_alias, '--', c='green', label='No ADAA')
+plt.plot(freqs_os, fft_os, '--', c='orange', label=f'OS{OS}')
+plt.plot(freqs_ad1, fft_ad1, 'green', label=f'ADAA1 + OS{OS}')
+plt.plot(freqs_ad2, fft_ad2, 'blue', label=f'ADAA2 + OS{OS}')
+plt.legend()
+plt.scatter(freqs_analog[peak_idxs], fft_analog[peak_idxs], c='r', marker='x')
+plt.xlim(0, 20000)
+plt.ylim(5)
+plt.title('Hard Clipping Distortion')
+plt.ylabel('Magnitude [dB]')
+plt.xlabel('Frequency [Hz]')
+plt.grid()
+
+plt.show()

scripts/anc.py

@@ -0,0 +1,122 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+
+def snr(output1, output2, reference):
+    return np.mean(reference ** 2) / np.mean((output1 - output2) ** 2)
+
+def normalize_2d(matrix):
+    norm = np.linalg.norm(matrix)
+    matrix = matrix/norm  # normalized matrix
+    return matrix
+
+
+def anc_adaptive_filter(audio):
+    # Initialize the adaptive filter
+    filter_length = 31
+    weights = np.zeros(filter_length)
+
+    # Run the adaptive filter
+    output = np.zeros_like(audio)
+    input_history = np.zeros_like(weights)
+    for i in range(len(audio)):
+        # Shift the input history and append the current input sample
+        input_history = np.roll(input_history, -1)
+        input_history[-1] = audio[i]
+
+        # Compute the filter output
+        output[i] = weights.dot(input_history)
+
+        # Update the filter weights
+        error = audio[i] - output[i]
+        weights += 0.1 * error * input_history
+
+    # Subtract the filtered noise from the original signal to obtain the noise-cancelled signal
+    # noise_cancelled = audio - output
+    return output
+
+
+def anc_kalman_filter(signal):
+    # Define the system state
+    state = np.zeros((2, 1))  # initial state (position and velocity)
+
+    # Define the system matrix
+    A = np.array([[1, 1], [0, 1]])  # state transition matrix
+
+    # Define the process noise covariance
+    Q = np.array([[0.001, 0.001], [0.001, 0.01]])  # covariance matrix
+
+    # Define the measurement matrix
+    H = np.array([[1, 0]])  # measurement matrix
+
+    # Define the measurement noise covariance
+    R = np.array([[0.1]])  # covariance matrix
+
+    # Define the Kalman gain
+    K = np.zeros((2, 1))  # initial Kalman gain
+
+    # Define the estimation error covariance
+    P = np.array([[1, 0], [0, 1]])  # initial estimation error covariance
+
+    # Initialize the output signal
+    output = np.zeros_like(signal)
+
+    # Iterate over the time steps
+    for i in range(1, len(signal)):
+        # Predict the state
+        state_pred = A @ state  # matrix multiplication
+        P_pred = A @ P @ A.T + Q  # covariance prediction
+
+        # Get the measurement
+        measurement = signal[i]
+
+        # Update the Kalman gain
+        K = P_pred @ H.T @ np.linalg.inv(H @ P_pred @ H.T + R)
+
+        # Update the state estimate
+        state = state_pred + K @ (measurement - H @ state_pred)
+
+        # Update the estimation error covariance
+        P = (np.eye(2) - K @ H) @ P_pred
+
+        # Store the output sample
+        output[i] = H @ state
+
+
+def gen_signal_with_noise(sample_rate, duration):
+    t = np.linspace(0, duration, int(duration * sample_rate), endpoint=False)
+    signal = np.sin(2 * np.pi * 440 * t)*0.5
+    signal += np.sin(2 * np.pi * 220 * t)*0.5
+    signal += np.sin(2 * np.pi * 110 * t)*0.25
+    noise = np.random.normal(scale=0.25, size=len(t))
+    return (signal, noise, t)
+
+
+# Generate a synthetic audio signal with noise
+sample_rate = 22050
+duration = 0.25
+
+signal, noise, t = gen_signal_with_noise(sample_rate, duration)
+noisy_signal = signal+noise*1.1
+output = anc_adaptive_filter(noisy_signal)
+
+# print(f'SNR: {snr(signal,noise,signal)}')
+print(f'SNR: {snr(noisy_signal,signal,noise)}')
+print(f'SNR: {snr(output,signal,noise)}')
+# Plot the input and output signals
+plt.plot(t, noisy_signal, label='Noisy')
+plt.plot(t, output, label='Clean')
+plt.xlabel('Time (s)')
+plt.ylabel('Amplitude')
+plt.legend()
+plt.show()
+
+# Plot the spectrum of the signal
+frequencies = np.fft.rfftfreq(len(noisy_signal), 1.0 / sample_rate)
+plt.semilogx(frequencies, np.abs(np.fft.rfft(normalize_2d(noisy_signal))) / sample_rate,label='Noisy')
+plt.semilogx(frequencies, np.abs(np.fft.rfft(normalize_2d(output))) / sample_rate,label='Clean')
+plt.gca().set_yscale('log')
+plt.xlabel("Normalized frequency")
+plt.ylabel("Magnitude")
+plt.legend()
+plt.show()

scripts/bandlimit.py

@@ -0,0 +1,78 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+
+def sawtooth(sample_rate, base_frequency):
+    # Calculate the period of the sawtooth wave
+    period = 1.0 / base_frequency
+
+    # Calculate the length of the wavetable in samples
+    wavetable_length = 2048
+
+    # Calculate the actual base frequency of the wavetable
+    wavetable_frequency = sample_rate / wavetable_length
+
+    # Generate a sawtooth wavetable with the desired length
+    t = np.linspace(0, period, wavetable_length, endpoint=False)
+    return t * wavetable_frequency - np.mod(t, 1.0 / wavetable_frequency)
+
+
+def sin_with_noise(sample_rate, frequency, noise_amplitude=0.1):
+    signal = np.sin(2 * np.pi * frequency * np.linspace(0, 1, sample_rate))
+    noise = noise_amplitude * np.random.randn(len(signal))
+    signal += noise
+    return signal
+
+
+def bandlimit(sample_rate, signal, low, high):
+    # Compute the discrete Fourier transform of the wavetable
+    dft = np.fft.rfft(signal)
+
+    # Normalize the frequencies
+    frequencies = np.fft.rfftfreq(len(wavetable), 1.0 / sample_rate)
+    normalized_frequencies = frequencies / (sample_rate / 2)
+
+    # Set the coefficients corresponding to frequencies outside the desired range to zero
+    low_frequency = low
+    high_frequency = high
+    dft[normalized_frequencies < low_frequency] = 0
+    dft[normalized_frequencies > high_frequency] = 0
+
+    # Compute the inverse discrete Fourier transform of the modified wavetable
+    filtered_wavetable = np.fft.irfft(dft)
+    return filtered_wavetable
+
+
+# Define the wavetable and sample rate
+sample_rate = 44100
+# wavetable = sawtooth(sample_rate, 440.0)
+wavetable = sin_with_noise(sample_rate, 440.0, 0.15)
+filtered_wavetable = bandlimit(sample_rate, wavetable, 0, 1.0)
+
+
+# Print the original and filtered wavetables
+print("Original wavetable:", wavetable)
+print("Filtered wavetable:", filtered_wavetable)
+
+# Plot the original and filtered wavetables
+plt.plot(wavetable, label="Original")
+plt.plot(filtered_wavetable, label="Filtered")
+
+plt.legend()
+plt.show()
+
+
+# # Normalize the frequencies
+frequencies = np.fft.rfftfreq(len(wavetable), 1.0 / sample_rate)
+# normalized_frequencies = frequencies / (sample_rate / 2)
+
+# Plot the spectrum of the wavetable
+plt.semilogx(frequencies[4:], np.abs(np.fft.rfft(wavetable))[4:]/sample_rate)
+plt.semilogx(frequencies[4:], np.abs(
+    np.fft.rfft(filtered_wavetable))[4:]/sample_rate)
+plt.xlabel("Normalized frequency")
+plt.ylabel("Magnitude")
+# Set the tick locations and labels on the x-axis
+ticks = [20, 55, 110, 220, 440, 880, 1760, 3520, 7040, 10000]
+plt.xticks(ticks, ticks)
+plt.show()

scripts/conv.py

@@ -0,0 +1,39 @@
+import numpy as np
+
+# Generate a synthetic audio signal with noise
+sample_rate = 44100
+duration = 1.0
+t = np.linspace(0, duration, int(duration * sample_rate), endpoint=False)
+signal = np.sin(2 * np.pi * 440 * t) + np.random.normal(scale=0.1, size=len(t))
+
+# Design a lowpass FIR filter
+cutoff_frequency = 1000  # Hz
+filter_length = 31
+window = 'hamming'
+fir_filter = signal.firwin(filter_length, cutoff_frequency, nyq=sample_rate / 2, window=window)
+
+# Set the frame size and overlap
+frame_size = 1024
+overlap = 512
+
+# Compute the number of frames
+num_frames = len(signal) // overlap - 1
+
+# Initialize the output signal
+output = np.zeros_like(signal)
+
+# Filter the signal using overlapping frames
+for i in range(num_frames):
+    # Get the current frame
+    start = i * overlap
+    stop = start + frame_size
+    frame = signal[start:stop]
+
+    # Apply the FIR filter to the frame using convolution
+    filtered_frame = np.convolve(frame, fir_filter, mode='valid')
+
+    # Overlap and add the filtered frame to the output signal
+    output[start:start+len(filtered_frame)] += filtered_frame
+
+# Normalize the output signal
+output /= np.max(np.abs(output))