Reconstructing Audio From Spectrogram (Using it as a filter)

Question

Given a spectrogram calculated using the following code:

import matplotlib.pyplot as plt
import numpy as np
import scipy
from scipy import signal, fft
from scipy.io import wavfile
from skimage import util
import librosa
import pylab
# from scipy.fftpack import fft

def stft(x, fs, framesz, hop):
    framesamp = int(framesz*fs)
    hopsamp = int(hop*fs)
    w = scipy.hanning(framesamp)
    X = scipy.array([scipy.fft(w*x[i:i+framesamp]) for i in range(0, len(x)-framesamp, hopsamp)])
    return X

def istft(X, fs, T, hop):
    x = scipy.zeros(T*fs)
    framesamp = X.shape[1]
    hopsamp = int(hop*fs)
    for n,i in enumerate(range(0, len(x)-framesamp, hopsamp)): x[i:i+framesamp] += scipy.real(scipy.ifft(X[n]))
    return x

audio_data = f'../recordings/SDRuno_20220309_140138Z_1kHz.wav' # y movement as test...
fs, data = wavfile.read(audio_data)
data = np.mean(data, axis=1) # convert to mono by avg l+r channels
N = data.shape[0]
L = N / fs
amp = 2 / N * np.abs(data)

print(f'Audio length: {L:.2f} seconds')

# f, ax = plt.subplots()
# ax.plot(np.arange(N) / fs, data)
# ax.set_xlabel('Time [s]')
# ax.set_ylabel('Amplitude [unknown]');
# plt.show()

M = int(fs * 0.001 * 20) # originally we had 1024, however now we use this for 20ms window resolution
# amp = 2 * np.sqrt(2)
stft_sig_f, stft_sig_t, stft_sig_Zxx = signal.stft(data, fs=fs, window='hann', nperseg=M, detrend=False)
# plt.pcolormesh(stft_sig_t, stft_sig_f, np.abs(stft_sig_Zxx), vmin=0, vmax=amp, shading='gouraud')
# plt.title('STFT Magnitude')
# plt.ylabel('Frequency [Hz]')
# plt.xlabel('Time [sec]')
# plt.show()

freqs, times, Sx = signal.spectrogram(data, fs=fs, window='hanning', nperseg=M, detrend=False, scaling='spectrum')
f, ax = plt.subplots(figsize=(10,5))
ax.pcolormesh(times, freqs / 1000, 10 * np.log10(Sx), cmap='inferno')
ax.set_ylabel('Frequency [kHz]')
ax.set_xlabel('Time [s]');
plt.show()

How can one use the spectrogram's "cleaner signal" to filter the audio. I did read somewhere that the Sx is effectively the Zxx**2 for which the Zxx can be used for an ISTFT to reconstruct audio, however my knowledge is pretty lacking. Any help would be appreciated!

EDIT: The current code for plotting the STFT just shows a black window.

The IQ WAV file used in this example can be found here: https://www.dropbox.com/s/at31myl5oufvfa3/SDRuno_20220309_140138Z_1kHz.wav?dl=0

Answer 1

"Spectrogram isn't invertible" isn't entirely true. Simply Google "invert spectrogram python" - among the results should be the most well-known algorithm, Griffin-Lim. One can hence tweak the spectrogram and then invert it, as a form of non-linear filtering. However, the best we can do is recover the signal within a global phase shift, $e^{j\phi_0}$ (which still preserves instantaneous frequency and amplitude perfectly).

Better results can be obtained by implementing STFT via differentiable operators and minimizing reconstruction error via gradient descent (scalogram example); torch.stft might help here.

Answer 2

The magnitude only contains half the required information, as @dan mentions.

You may try to set the complex phase angle to random values [0…2*pi]. That gives a «vocoder-ish» sound. For speech signals it sounds like «noise speaking».

In some cases phase is assumed to be minimum phase but I dont know if that makes sense in a spectrogram?

edit:

Perhaps I do not use the proper terminology. The below piece of MATLAB code performs an fft-based filterbank with 50% overlap, power-complimentary (dual) windowing and complex representation of fft coefficients. It is invertible (within numerical precision)

%% generate window
win_pow = 2;%power complementary
Nh = 256;
N = 2*Nh;
win = hann(N+1);
win = win(1:N);
win = win / sum(win([1 end/2+1]));
win = win.^(1/win_pow);
assert(max(abs(win.^win_pow + circshift(win, Nh).^win_pow - 1)) < 1e-14);

%% test 50% overlapping STFT
x = 2*rand(3*N,1)-1;
X = buffer(x, N, Nh, "nodelay");
X = X .* win;
X = fft(X);
% <- this is where you would do something useful
Y = ifft(X);
Y = Y .* win;
x_hat = debuffer(Y, Nh);
err = x(1:length(x_hat)) - x_hat;
assert(max(abs(err((Nh+1):2*N))) < 1e-15);

This filterbank is non-overlapping, rectangular windowed and also invertible. I thought that this, too was a "STFT" if somewhat poor due to being "critically sampled":

%% test critically sampled STFT
x = 2*rand(3*N,1)-1;
X = buffer(x, N, 0, "nodelay");
X = X;
X = fft(X);
Y = ifft(X);
Y = Y;
x_hat = debuffer(Y, 0);
err = x(1:length(x_hat)) - x_hat;
assert(max(abs(err((Nh+1):2*N))) < 1e-15);

Whenever I see abs(X) (and the phase is thrown away), I think of that as a spectrogram, and using the critically sampled filterbank from the previous snippet, I cannot see how it would be invertible without some significant assumptions about signal statistics? Is this what you guys call a STFT? A Spectrogram? I am confused by the responses.

%% test critically sampled spectrogram
x = 2*rand(3*N,1)-1;
X = buffer(x, N, 0, "nodelay");
X = X;
X = fft(X);
X = abs(X);
Y = ifft(X);
Y = Y;
x_hat = debuffer(Y, 0);
err = x(1:length(x_hat)) - x_hat;

x is represented as a 1536x1 vector, 12288 bytes on my machine.

X is represented as a 512x3 array, 12288 bytes on my machine. However, as it is mirrored, it only contains 257x3x8 = 6168 bytes of unique information. Thus it cannot possibly convey all possibilities of the x vector?

I see how an oversampled filterbank would contain "additional information" that (even if you throw away phase information) might be used to aid reconstruction. It sounds like an exercise I would rather not venture into.