This Python script provides an implementation of the Synchrosqueezing Transform (SST) for one-dimensional signals. The code is a direct translation of a Julia script, designed to preserve the specific phase relationships required for the algorithm to work correctly.
The Synchrosqueezing Transform is a time-frequency analysis technique that sharpens the time-frequency representation (TFR) of a signal, making it easier to identify and extract individual frequency components, especially for signals where the frequency changes over time.
The script requires the following Python libraries:
numpyscipymatplotlib(for the example usage)
You can install these dependencies using pip:
pip install numpy scipy matplotlibThis is the main function that computes the Synchrosqueezing Transform.
x(1D numpy array): The input signal.fs(int, optional): The sampling frequency of the signal. Defaults to 1.t(1D numpy array, optional): The time steps at which to compute the STFT. If not provided, it defaults tonp.arange(len(x)).window(str, optional): The window function to use. Defaults to 'hann'. Must be a window with an odd length.nperseg(int, optional): The number of samples per segment of the STFT. Defaults to 256 or the length of the signal, whichever is smaller.nfft(int, optional): The number of FFT points. Defaults to the length of the signal.
f(1D numpy array): The frequency bins.t(1D numpy array): The time steps.TFR(2D numpy array): The standard Short-Time Fourier Transform (STFT) of the signal.RTFR(2D numpy array): The Synchrosqueezed STFT of the signal.
The script can be run directly from the command line to see an example of the Synchrosqueezing Transform applied to a chirp signal.
python sq_stft.pyThis will generate a plot showing two subplots:
- The standard Short-Time Fourier Transform (STFT).
- The Synchrosqueezed STFT (sqSTFT).
The red dashed line in both plots indicates the true instantaneous frequency of the chirp signal. As you can see, the Synchrosqueezed STFT provides a much sharper and more concentrated representation of the signal's frequency content over time compared to the standard STFT.
import matplotlib.pyplot as plt
from scipy.signal import chirp
# 1. Generate a chirp signal
fs = 1000
T = 2
time_vec = np.linspace(0, T, int(fs * T), endpoint=False)
f_start = 10
f_end = fs / 4
sig = chirp(time_vec, f0=f_start, f1=f_end, t1=T, method='linear')
# 2. Compute STFT and Synchrosqueezed STFT
n_fft = 512
n_perseg = 129 # Must be odd
hop_size = 2
t_steps = np.arange(0, len(sig), hop_size)
f, times, TFR, RTFR = sqSTFT(sig, fs=fs, t=t_steps, window='hann', nperseg=n_perseg, nfft=n_fft)
# 3. Plot the results
# ... (plotting code) ...This Python script provides an implementation of the Inverse Synchrosqueezing Transform (i-sqSTFT) for one-dimensional signals. It works in conjunction with the sq_stft.py script, which provides the forward Synchrosqueezing Transform. The inverse transform aims to reconstruct the original time-domain signal from its Synchrosqueezed representation.
The reconstruction method implemented here leverages both the original Short-Time Fourier Transform (STFT) matrix (TFR) and the Synchrosqueezed STFT matrix (RTFR). The RTFR is used to create a mask that highlights the dominant frequency components. This mask is then applied to the TFR, preserving the crucial phase information from the original STFT only for the synchrosqueezed components. Finally, a standard Inverse Short-Time Fourier Transform (ISTFT) is performed on this masked TFR to recover the time-domain signal.
The script requires the following Python libraries:
numpyscipymatplotlib(for the example usage)sq_stft(the forward Synchrosqueezing Transform script, which must be in the same directory or accessible via Python path)
You can install numpy, scipy, and matplotlib using pip:
pip install numpy scipy matplotlibEnsure that sq_stft.py is present in the same directory as inverse_sqSTFT.py for successful execution.
This is the main function that computes the Inverse Synchrosqueezing Transform.
TFR(2D numpy array): The original STFT matrix (positive frequencies) obtained fromsqSTFT.RTFR(2D numpy array): The synchrosqueezed STFT matrix (positive frequencies) obtained fromsqSTFT.fs(float): The sampling frequency of the original signal.t(1D numpy array): The time steps used in the forward STFT.window(str): The name of the window function used in the forward STFT (e.g., 'hann').nperseg(int): The number of samples per segment used in the forward STFT.nfft(int): The number of points for the FFT used in the forward STFT.
np.ndarray: The reconstructed time-domain signal.
The script can be run directly from the command line to see an example of the inverse Synchrosqueezing Transform applied to a chirp signal.
python inverse_sqSTFT.pyThis will generate several plots:
- Original Signal: The synthetic chirp signal used for the demonstration.
- Reconstructed Signal (from i-sqSTFT): The signal recovered after applying the forward and inverse Synchrosqueezing Transforms.
- Difference (Error): The point-wise difference between the original and reconstructed signals, illustrating the accuracy of the reconstruction.
- Synchrosqueezed Representation (RTFR): The time-frequency representation after the forward Synchrosqueezing Transform.
- Masked Original STFT (TFR): The original STFT with the mask derived from the
RTFRapplied, showing which components are used for reconstruction.
These plots collectively demonstrate the process and the fidelity of the signal reconstruction using the inverse Synchrosqueezing Transform.
import numpy as np
from scipy.signal import chirp
import matplotlib.pyplot as plt
from sq_stft import sqSTFT # Assuming sq_stft.py is available
# 1. Create a synthetic signal
fs = 1000
T = 2
time_vec = np.linspace(0, T, int(fs * T), endpoint=False)
f_start = 10
f_end = fs / 4
sig = chirp(time_vec, f0=f_start, f1=f_end, t1=T, method='linear')
sig = sig / np.max(np.abs(sig))
# 2. Define parameters for transforms
n_fft = 512
n_perseg = 129
hop_size = 2
t_steps = np.arange(0, len(sig), hop_size)
window_type = 'hann'
# 3. Perform forward sqSTFT
f, times, TFR, RTFR = sqSTFT(sig, fs=fs, t=t_steps, window=window_type, nperseg=n_perseg, nfft=n_fft)
# 4. Perform inverse sqSTFT
reconstructed_signal = i_sqSTFT(TFR, RTFR, fs=fs, t=t_steps, window=window_type, nperseg=n_perseg, nfft=n_fft)
# ... (plotting and comparison code) ...