From a5f7b156627ce22a05e78c3fd97ac6a85c4de08e Mon Sep 17 00:00:00 2001 From: Tobias Hienzsch Date: Mon, 5 Feb 2024 21:15:48 +0100 Subject: [PATCH] Copy code examples from mc-audio-effects --- requirements.txt | 2 + scripts/adaa.py | 133 +++++++++++++++++++++++++++++++++++++ scripts/anc.py | 122 ++++++++++++++++++++++++++++++++++ scripts/bandlimit.py | 78 ++++++++++++++++++++++ scripts/conv.py | 39 +++++++++++ scripts/dk_method.py | 35 ++++++++++ scripts/dk_method_diode.py | 32 +++++++++ scripts/integral.py | 77 +++++++++++++++++++++ scripts/noise.py | 38 +++++++++++ scripts/phase_vocoder.py | 68 +++++++++++++++++++ scripts/stretch.py | 28 ++++++++ scripts/stretch_test.py | 117 ++++++++++++++++++++++++++++++++ scripts/svf_peak.cpp | 80 ++++++++++++++++++++++ scripts/svf_response.py | 35 ++++++++++ 14 files changed, 884 insertions(+) create mode 100644 scripts/adaa.py create mode 100644 scripts/anc.py create mode 100644 scripts/bandlimit.py create mode 100644 scripts/conv.py create mode 100644 scripts/dk_method.py create mode 100644 scripts/dk_method_diode.py create mode 100644 scripts/integral.py create mode 100644 scripts/noise.py create mode 100644 scripts/phase_vocoder.py create mode 100644 scripts/stretch.py create mode 100644 scripts/stretch_test.py create mode 100644 scripts/svf_peak.cpp create mode 100644 scripts/svf_response.py diff --git a/requirements.txt b/requirements.txt index f351582..ecc379f 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,4 +1,6 @@ jupyter matplotlib numpy +scipy +sympy tabulate diff --git a/scripts/adaa.py b/scripts/adaa.py new file mode 100644 index 0000000..797e82e --- /dev/null +++ b/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() diff --git a/scripts/anc.py b/scripts/anc.py new file mode 100644 index 0000000..d0ac7b2 --- /dev/null +++ b/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() diff --git a/scripts/bandlimit.py b/scripts/bandlimit.py new file mode 100644 index 0000000..66e8539 --- /dev/null +++ b/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() diff --git a/scripts/conv.py b/scripts/conv.py new file mode 100644 index 0000000..93f7e3b --- /dev/null +++ b/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)) \ No newline at end of file diff --git a/scripts/dk_method.py b/scripts/dk_method.py new file mode 100644 index 0000000..a51c1e3 --- /dev/null +++ b/scripts/dk_method.py @@ -0,0 +1,35 @@ +import numpy as np + +# Circuit parameters +R = 10.0 # Resistance (Ohms) +C = 0.01 # Capacitance (Farads) + +# Nonlinear component parameters +V_threshold = 1.0 # Threshold voltage (Volts) +I_max = 1.0 # Maximum current (Amps) + +# Initial conditions +v_in = 5.0 # Input voltage (Volts) +v_out = 0.0 # Output voltage (Volts) +i_c = 0.0 # Capacitor current (Amps) + +# Time step +dt = 0.001 # Time step (Seconds) + +# Number of time steps +num_steps = 100 + +# Compute the state of the circuit at each time step +for t in range(num_steps): + # Compute the capacitor current + i_c += dt * (v_in - v_out) / (R * C) + # Compute the output voltage + v_out = v_in + i_c * R + # Compute the current through the nonlinear component + i_nl = np.clip(i_c, 0.0, I_max) + # Compute the voltage drop across the nonlinear component + v_nl = V_threshold * (1.0 - np.exp(-i_nl)) + # Compute the final output voltage + v_out += v_nl + # Print the output voltage + print(v_out) \ No newline at end of file diff --git a/scripts/dk_method_diode.py b/scripts/dk_method_diode.py new file mode 100644 index 0000000..019c6d4 --- /dev/null +++ b/scripts/dk_method_diode.py @@ -0,0 +1,32 @@ +import numpy as np + +# Circuit parameters +R = 100.0 # Resistance of resistor (Ohms) +C = 0.01 # Capacitance of capacitor (Farads) +V_T = 0.0257 # Thermal voltage (Volts) + +# Initial conditions +v_in = 5.2 # Input voltage (Volts) +v_out = 5.0 # Output voltage (Volts) +i_c = 0.0 # Capacitor current (Amps) + +# Time step +dt = 0.001 # Time step (Seconds) + +# Number of time steps +num_steps = 10 + +# Compute the state of the circuit at each time step +for t in range(num_steps): + # Compute the capacitor current + i_c += dt * (v_in - v_out) / (R * C) + # Compute the output voltage + v_out = v_in + i_c * R + # Compute the diode current + i_d = np.clip(i_c, 0.0, np.inf) + # Compute the voltage drop across the diode + v_d = V_T * np.log(1.0 + i_d / V_T) + # Compute the final output voltage + v_out += v_d + # Print the output voltage + print(v_out) diff --git a/scripts/integral.py b/scripts/integral.py new file mode 100644 index 0000000..5edfed5 --- /dev/null +++ b/scripts/integral.py @@ -0,0 +1,77 @@ +import matplotlib.pyplot as plt +import numpy as np +import sympy +from sympy.utilities.lambdify import lambdify + +sympy.init_printing(use_unicode=False, wrap_line=False) + +# x = sympy.Symbol('x', positive=True, nonzero=True) +x = sympy.Symbol('x', real=True) + + +tanh_clip = sympy.tanh(x) +sin_clip = -sympy.sin(sympy.pi*0.5*x) + + +degree = 3 +one_over_degree = sympy.Rational(1, degree) +norm_factor = sympy.Rational(degree - 1, degree) +inv_norm = sympy.Rational(1, norm_factor) +y = x * norm_factor +soft3_clip = (y - (y**degree) * one_over_degree)*inv_norm + +degree = 5 +one_over_degree = sympy.Rational(1, degree) +norm_factor = sympy.Rational(degree - 1, degree) +inv_norm = sympy.Rational(1, norm_factor) +y = x * norm_factor +soft5_clip = (y - (y**degree) * one_over_degree)*inv_norm + +alpha = sympy.Symbol('alpha') +beta = sympy.Symbol('beta') +alpha = 1.79 # sympy.Rational(179, 100) +beta = sympy.Rational(1, 5) +diode_clip = beta * (sympy.exp(alpha*x) - 1.0) + + +# x_in = np.linspace(-1.25, 1.25, 128) +# sin_out = lambdify(x, sin_clip, 'numpy')(x_in) +# tanh_out = lambdify(x, tanh_clip, 'numpy')(x_in) +# soft3_out = lambdify(x, soft3_clip, 'numpy')(x_in) +# soft5_out = lambdify(x, soft5_clip, 'numpy')(x_in) +# diode_out = lambdify(x, diode_clip, 'numpy')(x_in) + + +# plt.plot(x_in, tanh_out, label="tanh") +# plt.plot(x_in, sin_out, label="sin") +# plt.plot(x_in, soft3_out, label="soft3") +# plt.plot(x_in, soft5_out, label="soft5") +# plt.plot(x_in, diode_out, label="diode") +# plt.grid(which="major", linewidth=1) +# plt.grid(which="minor", linewidth=0.2) +# plt.minorticks_on() +# plt.legend() +# plt.show() + + +f = sympy.sign(x)*2*x +f = (3-(2-3*x)**2) / 3 +f = 1 +f = sympy.Piecewise( # Hard-clip + (-1 , sympy.Le(x, -1)), + (1 , sympy.Ge(x, 1)), + (x, True), +) +# f = sympy.Piecewise( # Overdrive +# (sympy.sign(x)*2*x, sympy.And(sympy.Ge(x, 0) , sympy.Le(x, sympy.Rational(1,3)))), +# ((3-(2-3*x)**2) / 3, sympy.Gt(x, sympy.Rational(1,3)) & sympy.Le(x, sympy.Rational(2,3))), +# (1, True) +# ) +# f = sympy.sign(x)*(1-sympy.exp(-sympy.Abs(x))) +# f = sympy.Abs(x) # Full-wave +# f = sympy.Piecewise((x, sympy.Gt(x, 0)), (sympy.Abs(x), True)) # Half-wave +AD1 = sympy.integrate(f, x) +AD2 = sympy.integrate(AD1, x) +print(f) +print(AD1) +print(AD2) diff --git a/scripts/noise.py b/scripts/noise.py new file mode 100644 index 0000000..be5688a --- /dev/null +++ b/scripts/noise.py @@ -0,0 +1,38 @@ +import numpy as np + + +def test_roundtrip_fft_ift(): + # Generate a test signal + audio = np.sin(2 * np.pi * 440.0 * np.linspace(0, 1, 44100)) + + # Add white Gaussian noise to the signal + noise_amplitude = 0.1 + noise = noise_amplitude * np.random.randn(len(audio)) + audio += noise + + # Compute the FFT of the test signal + fft = np.fft.fft(audio) + + # Compute the IFFT of the FFT + ifft = np.fft.ifft(fft) + + # Compute the relative error between the original and reconstructed signals + relative_error = np.abs(audio - ifft) / np.abs(audio) + + # Compute the RMS error between the original and reconstructed signals + rms_error = np.sqrt(np.mean((audio - ifft) ** 2)) + + # Compute the PSNR between the original and reconstructed signals + psnr = 10 * np.log10(np.max(audio) ** 2 / np.mean((audio - ifft) ** 2)) + + # Compute the SNR between the original and reconstructed signals + snr = np.mean(audio ** 2) / np.mean((audio - ifft) ** 2) + + # Print the results + # print("Relative error:", relative_error) + print("RMS error:", np.abs(rms_error)) + print("PSNR:", psnr) + print("SNR:", snr) + + +test_roundtrip_fft_ift() diff --git a/scripts/phase_vocoder.py b/scripts/phase_vocoder.py new file mode 100644 index 0000000..e0121cd --- /dev/null +++ b/scripts/phase_vocoder.py @@ -0,0 +1,68 @@ +import numpy as np + + +def phase_vocoder(audio, rate, time_stretch_ratio, block_size=1024, window_type='hann'): + # Create the window function + window = create_window(block_size, window_type) + + # Calculate the time-stretch factor + stretch_factor = rate / (rate * time_stretch_ratio) + + # Initialize the output audio + stretched_audio = np.zeros(len(audio)) + + # Process the audio in blocks + for i in range(0, len(audio), block_size): + # Get the current block of audio + block = audio[i:i+block_size] + + # Apply the window function to the block + windowed_block = block * window + + # Convert the block to the frequency domain + spectrum = np.fft.rfft(windowed_block) + + # Get the length of the spectrum + M = len(spectrum) + + # Create an empty array to store the stretched spectrum + stretched_spectrum = np.zeros(M, dtype=complex) + + # Stretch the spectrum by interpolating the magnitude and phase + for j in range(M): + # Calculate the stretched frequency index + k = int(j * stretch_factor) + + # Interpolate the magnitude and phase + if k < M: + stretched_spectrum[j] = spectrum[k] + else: + stretched_spectrum[j] = 0 + + # Convert the stretched spectrum back to the time domain + stretched_block = np.fft.irfft(stretched_spectrum) + + # Overlap-add the stretched block to the output audio + stretched_audio[i:i+block_size] += stretched_block + + # Normalize the audio + stretched_audio /= np.max(np.abs(stretched_audio)) + + return stretched_audio + + +def create_window(N, window_type='hann'): + # Create a Hann window + if window_type == 'hann': + window = 0.5 * (1 - np.cos(2 * np.pi * np.arange(N) / (N - 1))) + # Create a Hamming window + elif window_type == 'hamming': + window = 0.54 - 0.46 * np.cos(2 * np.pi * np.arange(N) / (N - 1)) + # Create a Blackman window + elif window_type == 'blackman': + window = 0.42 - 0.5 * np.cos(2 * np.pi * np.arange(N) / (N - 1)) + \ + 0.08 * np.cos(4 * np.pi * np.arange(N) / (N - 1)) + else: + raise ValueError('Invalid window type') + + return window diff --git a/scripts/stretch.py b/scripts/stretch.py new file mode 100644 index 0000000..f5c6fd4 --- /dev/null +++ b/scripts/stretch.py @@ -0,0 +1,28 @@ +import numpy as np + + +def stretch(audio, stretch_factor, block_size, window): + # Initialize the output array + stretched_audio = np.zeros_like(audio) + + # Iterate over blocks of the audio signal + for i in range(0, len(audio), block_size): + # Window the current block + block = audio[i:i+block_size] * window + + # Compute the STFT of the current block + stft = np.fft.fft(block) + + # Shift the STFT to center the frequency spectrum + stft = np.fft.fftshift(stft) + + # Compute the stretched STFT + start = int(len(stft) * (1 - stretch_factor) / 2) + stop = int(len(stft) * (1 + stretch_factor) / 2) + stretched_stft = np.zeros_like(stft) + stretched_stft[start:stop] = stft[start:stop] + + # Shift the stretched STFT back to its original position + stretched_stft = np.fft.ifftshift(stretched_stft) + + # Compute the inverse STFT of the stretched STFT diff --git a/scripts/stretch_test.py b/scripts/stretch_test.py new file mode 100644 index 0000000..fca4203 --- /dev/null +++ b/scripts/stretch_test.py @@ -0,0 +1,117 @@ +import pytest +import numpy as np +# import numpy as np + + +def stretch(audio, stretch_factor, block_size, window): + # Pad the audio signal with zeros if its length is not divisible by the block size + padding = block_size - len(audio) % block_size + audio = np.pad(audio, (0, padding), 'constant') + + # Initialize the output array + stretched_audio = np.zeros_like(audio) + + # Iterate over blocks of the audio signal + for i in range(0, len(audio), block_size): + # Window the current block + block = audio[i:i+block_size] * window + + # Compute the STFT of the current block + stft = np.fft.fft(block) + + # Shift the STFT to center the frequency spectrum + stft = np.fft.fftshift(stft) + + # Compute the stretched STFT + start = int(len(stft) * (1 - stretch_factor) / 2) + stop = int(len(stft) * (1 + stretch_factor) / 2) + stretched_stft = np.zeros_like(stft) + stretched_stft[start:stop] = stft[start:stop] + + # Shift the stretched STFT back to its original position + stretched_stft = np.fft.ifftshift(stretched_stft) + + # Compute the inverse STFT of the stretched STFT + stretched_block = np.fft.ifft(stretched_stft) + + # Overlap-add the stretched block onto the output array + stretched_audio[i:i+block_size] += np.real(stretched_block) + + # Trim the padded zeros from the output + stretched_audio = stretched_audio[:len(audio)-padding] + + # Return the stretched audio + return stretched_audio + +# from stretch_function import stretch + +# from .stretch import stretch + + +def test_stretch(): + # Generate a test signal + audio = np.sin(2 * np.pi * 440.0 * np.linspace(0, 1, 44100)) + + # Stretch the audio signal by a factor of 2 + stretch_factor = 2.0 + stretched_audio = stretch(audio, stretch_factor, 2048, np.hanning(2048)) + + # Check that the length of the stretched audio is correct + print(len(stretched_audio)) + print(len(audio) * stretch_factor) + assert len(stretched_audio) == len(audio) * stretch_factor + + +def test_stretch_minimum_block_size(): + # Generate a test signal + audio = np.sin(2 * np.pi * 440.0 * np.linspace(0, 1, 44100)) + + # Stretch the audio signal by a factor of 2 + stretch_factor = 2.0 + stretched_audio = stretch(audio, stretch_factor, 1, np.hanning(1)) + + # Check that the length of the stretched audio is correct + assert len(stretched_audio) == len(audio) * stretch_factor + + +def test_stretch_maximum_block_size(): + # Generate a test signal + audio = np.sin(2 * np.pi * 440.0 * np.linspace(0, 1, 44100)) + + # Stretch the audio signal by a factor of 2 + stretch_factor = 2.0 + stretched_audio = stretch(audio, stretch_factor, + len(audio), np.hanning(len(audio))) + + # Check that the length of the stretched audio is correct + assert len(stretched_audio) == len(audio) * stretch_factor + + +def test_stretch_zero_stretch_factor(): + # Generate a test signal + audio = np.sin(2 * np.pi * 440.0 * np.linspace(0, 1, 44100)) + + # Stretch the audio signal by a factor of 0 + stretch_factor = 0.0 + stretched_audio = stretch(audio, stretch_factor, 2048, np.hanning(2048)) + + # Check that the length of the stretched audio is correct + assert len(stretched_audio) == len(audio) * stretch_factor + + +def test_stretch_invalid_stretch_factor(): + # Generate a test signal + audio = np.sin(2 * np.pi * 440.0 * np.linspace(0, 1, 44100)) + + # Stretch the audio signal by a negative factor + stretch_factor = -1.0 + with pytest.raises(ValueError): + stretched_audio = stretch( + audio, stretch_factor, 2048, np.hanning(2048)) + + +test_stretch() +test_stretch_minimum_block_size() +test_stretch_maximum_block_size() +test_stretch_zero_stretch_factor() +test_stretch_invalid_stretch_factor() diff --git a/scripts/svf_peak.cpp b/scripts/svf_peak.cpp new file mode 100644 index 0000000..de82793 --- /dev/null +++ b/scripts/svf_peak.cpp @@ -0,0 +1,80 @@ +#include +#include +#include +#include + +static constexpr auto Pi = 3.14159265; + +inline float linearToDecibel(float amplitude) +{ + return 20.0f * std::log10(amplitude); +} + +inline float decibelToLinear(float dB) +{ + return std::pow(10.0f, dB / 20.0f); +} + +typedef struct +{ + double f0; // Center frequency of the peak + double q; // Quality factor of the peak + double fs; // Sample rate + double a, b, c; // Filter coefficients + double x1, x2; // Input state variables + double y1, y2; // Output state variables +} peak_filter_t; + +auto peak_filter_init(peak_filter_t *filter, double f0, double q, double fs) -> void +{ + filter->f0 = f0; + filter->q = q; + filter->fs = fs; + filter->x1 = 0; + filter->x2 = 0; + filter->y1 = 0; + filter->y2 = 0; +} + +auto peak_filter_process(peak_filter_t *filter, double *input, double *output, + uint32_t num_samples) -> void +{ + double omega = 2 * Pi * filter->f0 / filter->fs; + double alpha = sin(omega) / (2 * filter->q); + double a = 1 + alpha; + filter->a = alpha / a; + filter->b = 0; + filter->c = -alpha / a; + for (uint32_t i = 0; i < num_samples; ++i) + { + double x = input[i]; + double y = filter->a * x + filter->b * filter->x1 + filter->c * filter->x2 - + filter->b * filter->y1 - filter->c * filter->y2; + output[i] = y; + filter->x2 = filter->x1; + filter->x1 = x; + filter->y2 = filter->y1; + filter->y1 = y; + } +} + +auto peak_filter_response(peak_filter_t *filter, double f) -> std::complex +{ + using namespace std::complex_literals; + double omega = 2 * Pi * f / filter->fs; + std::complex z = exp(-1i * omega); + return (1.0 + filter->a * z + filter->c * z * z) / + (1.0 - filter->a * z + filter->c * z * z); +} + +auto main() -> int +{ + auto filter = peak_filter_t{}; + peak_filter_init(&filter, 440.0, 0.71, 44'100.0); + peak_filter_process(&filter, nullptr, nullptr, 0); + for (auto f : {20.0, 100.0, 200.0, 400.0, 440.0, 800.0, 1000.0, 5000.0}) + { + auto const z = peak_filter_response(&filter, f); + std::printf("f: %f, g: %f\n", f, linearToDecibel(std::abs(z))); + } +} \ No newline at end of file diff --git a/scripts/svf_response.py b/scripts/svf_response.py new file mode 100644 index 0000000..63af020 --- /dev/null +++ b/scripts/svf_response.py @@ -0,0 +1,35 @@ +import matplotlib.pyplot as plt +import numpy as np + +f = np.linspace(0.0, np.pi, 256) +z = np.exp(f*1.0j) +H01 = (z*0.1)/(z-(1-0.1)) +H02 = (z*0.2)/(z-(1-0.2)) +H05 = (z*0.5)/(z-(1-0.5)) +H10 = (z*1.0)/(z-(1-1.0)) +H11 = (z*1.1)/(z-(1-1.1)) + +plt.semilogx(f/(np.pi*2), 20*np.log10(np.abs(H01)), label="ω = 0.1") +plt.semilogx(f/(np.pi*2), 20*np.log10(np.abs(H02)), label="ω = 0.2") +plt.semilogx(f/(np.pi*2), 20*np.log10(np.abs(H05)), label="ω = 0.5") +plt.semilogx(f/(np.pi*2), 20*np.log10(np.abs(H10)), label="ω = 1.0") +plt.semilogx(f/(np.pi*2), 20*np.log10(np.abs(H11)), label="ω = 1.1") + +# plt.semilogx(f/(np.pi*2), np.angle(H01), label="ω = 0.1") +# plt.semilogx(f/(np.pi*2), np.angle(H02), label="ω = 0.2") +# plt.semilogx(f/(np.pi*2), np.angle(H05), label="ω = 0.5") +# plt.semilogx(f/(np.pi*2), np.angle(H10), label="ω = 1.0") +# plt.semilogx(f/(np.pi*2), np.angle(H11), label="ω = 1.1") + +# plt.vlines(0.02*np.pi,-25,0) +# plt.vlines(0.1*np.pi,-25,0) +# plt.vlines(np.pi,-25,0) + +plt.legend() +plt.grid(which="major", linewidth=1) +plt.grid(which="minor", linewidth=0.2) +plt.xticks([0.1, 0.5], [0.1, 0.5]) +plt.minorticks_on() +plt.show() + +# print(np.abs(H))