"""
Module: libfmp.c2.c2_fourier
Author: Frank Zalkow, Meinard Müller
License: The MIT license, https://opensource.org/licenses/MIT
This file is part of the FMP Notebooks (https://www.audiolabs-erlangen.de/FMP)
"""
import numpy as np
from numba import jit
import librosa
[docs]@jit(nopython=True)
def generate_matrix_dft(N, K):
"""Generates a DFT (discrete Fourier transfrom) matrix
Notebook: C2/C2_DFT-FFT.ipynb
Args:
N (int): Number of samples
K (int): Number of frequency bins
Returns:
dft (np.ndarray): The DFT matrix
"""
dft = np.zeros((K, N), dtype=np.complex128)
for n in range(N):
for k in range(K):
dft[k, n] = np.exp(-2j * np.pi * k * n / N)
return dft
[docs]@jit(nopython=True)
def generate_matrix_dft_inv(N, K):
"""Generates an IDFT (inverse discrete Fourier transfrom) matrix
Notebook: C2/C2_STFT-Inverse.ipynb
Args:
N (int): Number of samples
K (int): Number of frequency bins
Returns:
dft (np.ndarray): The IDFT matrix
"""
dft = np.zeros((K, N), dtype=np.complex128)
for n in range(N):
for k in range(K):
dft[k, n] = np.exp(2j * np.pi * k * n / N) / N
return dft
[docs]@jit(nopython=True)
def dft(x):
"""Compute the disrcete Fourier transfrom (DFT)
Notebook: C2/C2_DFT-FFT.ipynb
Args:
x (np.ndarray): Signal to be transformed
Returns:
X (np.ndarray): Fourier transform of x
"""
x = x.astype(np.complex128)
N = len(x)
dft_mat = generate_matrix_dft(N, N)/N
return np.dot(dft_mat, x)
[docs]@jit(nopython=True)
def idft(X):
"""Compute the inverse discrete Fourier transfrom (IDFT)
Args:
X (np.ndarray): Signal to be transformed
Returns:
x (np.ndarray): Inverse Fourier transform of X
"""
X = X.astype(np.complex128)
N = len(X)
dft_mat = generate_matrix_dft_inv(N, N)
return np.dot(dft_mat, X)
[docs]@jit(nopython=True)
def twiddle(N):
"""Generate the twiddle factors used in the computation of the fast Fourier transform (FFT)
Notebook: C2/C2_DFT-FFT.ipynb
Args:
N (int): Number of samples
Returns:
sigma (np.ndarray): The twiddle factors
"""
k = np.arange(N // 2)
sigma = np.exp(-2j * np.pi * k / N)
return sigma
[docs]@jit(nopython=True)
def twiddle_inv(N):
"""Generate the twiddle factors used in the computation of the Inverse fast Fourier transform (IFFT)
Args:
N (int): Number of samples
Returns:
sigma (np.ndarray): The twiddle factors
"""
n = np.arange(N // 2)
sigma = np.exp(2j * np.pi * n / N)
return sigma
[docs]@jit(nopython=True)
def fft(x):
"""Compute the fast Fourier transform (FFT)
Notebook: C2/C2_DFT-FFT.ipynb
Args:
x (np.ndarray): Signal to be transformed
Returns:
X (np.ndarray): Fourier transform of x
"""
x = x.astype(np.complex128)
N = len(x)
log2N = np.log2(N)
assert log2N == int(log2N), 'N must be a power of two!'
X = np.zeros(N, dtype=np.complex128)
if N == 1:
return x
else:
this_range = np.arange(N)
A = fft(x[this_range % 2 == 0])
B = fft(x[this_range % 2 == 1])
C = twiddle(N) * B
X[:N//2] = A + C
X[N//2:] = A - C
return X
[docs]@jit(nopython=True)
def ifft_noscale(X):
"""Compute the inverse fast Fourier transform (IFFT) without the final scaling factor of 1/N
Args:
X (np.ndarray): Fourier transform of x
Returns:
x (np.ndarray): Inverse Fourier transform of X
"""
X = X.astype(np.complex128)
N = len(X)
log2N = np.log2(N)
assert log2N == int(log2N), 'N must be a power of two!'
x = np.zeros(N, dtype=np.complex128)
if N == 1:
return X
else:
this_range = np.arange(N)
A = ifft_noscale(X[this_range % 2 == 0])
B = ifft_noscale(X[this_range % 2 == 1])
C = twiddle_inv(N) * B
x[:N//2] = A + C
x[N//2:] = A - C
return x
[docs]@jit(nopython=True)
def ifft(X):
"""Compute the inverse fast Fourier transform (IFFT)
Args:
X (np.ndarray): Fourier transform of x
Returns:
x (np.ndarray): Inverse Fourier transform of X
"""
return ifft_noscale(X) / len(X)
[docs]def stft_basic(x, w, H=8, only_positive_frequencies=False):
"""Compute a basic version of the discrete short-time Fourier transform (STFT)
Notebook: C2/C2_STFT-Basic.ipynb
Args:
x (np.ndarray): Signal to be transformed
w (np.ndarray): Window function
H (int): Hopsize (Default value = 8)
only_positive_frequencies (bool): Return only positive frequency part of spectrum (non-invertible)
(Default value = False)
Returns:
X (np.ndarray): The discrete short-time Fourier transform
"""
N = len(w)
L = len(x)
M = np.floor((L - N) / H).astype(int) + 1
X = np.zeros((N, M), dtype='complex')
for m in range(M):
x_win = x[m * H:m * H + N] * w
X_win = np.fft.fft(x_win)
X[:, m] = X_win
if only_positive_frequencies:
K = 1 + N // 2
X = X[0:K, :]
return X
[docs]def istft_basic(X, w, H, L):
"""Compute the inverse of the basic discrete short-time Fourier transform (ISTFT)
Notebook: C2/C2_STFT-Inverse.ipynb
Args:
X (np.ndarray): The discrete short-time Fourier transform
w (np.ndarray): Window function
H (int): Hopsize
L (int): Length of time signal
Returns:
x (np.ndarray): Time signal
"""
N = len(w)
M = X.shape[1]
x_win_sum = np.zeros(L)
w_sum = np.zeros(L)
for m in range(M):
x_win = np.fft.ifft(X[:, m])
# Avoid imaginary values (due to floating point arithmetic)
x_win = np.real(x_win)
x_win_sum[m * H:m * H + N] = x_win_sum[m * H:m * H + N] + x_win
w_shifted = np.zeros(L)
w_shifted[m * H:m * H + N] = w
w_sum = w_sum + w_shifted
# Avoid division by zero
w_sum[w_sum == 0] = np.finfo(np.float32).eps
x_rec = x_win_sum / w_sum
return x_rec, x_win_sum, w_sum
[docs]@jit(nopython=True)
def stft(x, w, H=512, zero_padding=0, only_positive_frequencies=False):
"""Compute the discrete short-time Fourier transform (STFT)
Args:
x (np.ndarray): Signal to be transformed
w (np.ndarray): Window function
H (int): Hopsize (Default value = 512)
zero_padding (bool): Number of zeros to be padded after windowing and before the Fourier transform of a frame
(Note: The purpose of this step is to increase the frequency sampling.) (Default value = 0)
only_positive_frequencies (bool): Return only positive frequency part of spectrum (non-invertible)
(Default value = False)
Returns:
X (np.ndarray): The discrete short-time Fourier transform
"""
N = len(w)
x = np.concatenate((np.zeros(N // 2), x, np.zeros(N // 2)))
L = len(x)
M = int(np.floor((L - N) / H)) + 1
X = np.zeros((N + zero_padding, M), dtype=np.complex128)
zero_padding_vector = np.zeros((zero_padding, ), dtype=x.dtype)
for m in range(M):
x_win = x[m * H:m * H + N] * w
if zero_padding > 0:
x_win = np.concatenate((x_win, zero_padding_vector))
X_win = fft(x_win)
# Note: X_win = np.fft.fft(x_win) does not work in combination with @jit
X[:, m] = X_win
if only_positive_frequencies:
K = 1 + (N + zero_padding) // 2
X = X[0:K, :]
return X
[docs]@jit(nopython=True)
def istft(X, w, H, L, zero_padding=0):
"""Compute the inverse discrete short-time Fourier transform (ISTFT)
Args:
X (np.ndarray): The discrete short-time Fourier transform
w (np.ndarray): Window function
H (int): Hopsize
L (int): Length of time signal
zero_padding (bool): Number of zeros to be padded after windowing and before the Fourier transform of a frame
(Default value = 0)
Returns:
x (np.ndarray): Reconstructed time signal
"""
N = len(w)
L = L + N
M = X.shape[1]
w_sum = np.zeros(L)
x_win_sum = np.zeros(L)
w_sum = np.zeros(L)
for m in range(M):
start_idx, end_idx = m * H, m * H + N + zero_padding
if start_idx > L:
break
x_win = ifft(X[:, m])
# Note: x_win = np.fft.ifft(X[:, m]) does not work in combination with @jit
if end_idx > L:
end_idx = L
x_win = x_win[:end_idx-start_idx]
cur_w = w[:end_idx-start_idx]
else:
cur_w = w
# Avoid imaginary values (due to floating point arithmetic)
x_win_real = np.real(x_win)
x_win_sum[start_idx:end_idx] = x_win_sum[start_idx:end_idx] + x_win_real
w_shifted = np.zeros(L)
w_shifted[start_idx:start_idx + len(cur_w)] = cur_w
w_sum = w_sum + w_shifted
# Avoid division by zero
w_sum[w_sum == 0] = np.finfo(np.float32).eps
x_rec = x_win_sum / w_sum
x_rec = x_rec[N // 2:-N // 2]
return x_rec
[docs]def stft_convention_fmp(x, Fs, N, H, pad_mode='constant', center=True, mag=False, gamma=0):
"""Compute the discrete short-time Fourier transform (STFT)
Notebook: C2/C2_STFT-FreqGridInterpol.ipynb
Args:
x (np.ndarray): Signal to be transformed
Fs (scalar): Sampling rate
N (int): Window size
H (int): Hopsize
pad_mode (str): Padding strategy is used in librosa (Default value = 'constant')
center (bool): Centric view as used in librosa (Default value = True)
mag (bool): Computes magnitude STFT if mag==True (Default value = False)
gamma (float): Constant for logarithmic compression (only applied when mag==True) (Default value = 0)
Returns:
X (np.ndarray): Discrete (magnitude) short-time Fourier transform
"""
X = librosa.stft(x, n_fft=N, hop_length=H, win_length=N,
window='hann', pad_mode=pad_mode, center=center)
if mag:
X = np.abs(X)**2
if gamma > 0:
X = np.log(1 + gamma * X)
F_coef = librosa.fft_frequencies(sr=Fs, n_fft=N)
T_coef = librosa.frames_to_time(np.arange(X.shape[1]), sr=Fs, hop_length=H)
# T_coef = np.arange(X.shape[1]) * H/Fs
# F_coef = np.arange(N//2+1) * Fs/N
return X, T_coef, F_coef