Matlab转Python：Stockwell变换函数ifft结果不一致问题排查

问题背景 #

我正在将Matlab的Stockwell变换函数stran转换为Python，测试数据保存为z.dat。Matlab版本的结果可视化正常，但Python代码生成的结果扭曲，与Matlab输出不一致，代码中除ifft部分外其余逻辑看似正常，需要排查问题。

Matlab绘图代码：

figure;
imagesc(abs(ST));
axis xy; 
xlabel("Time sample");
ylabel("Frequency");
title("S-transform result using stran.m")

注意：向stran.m传入的数据需为行向量。

用户的Python转换代码：

import numpy as np
from scipy import fftpack, linalg
from numpy import matlib
import matplotlib.pyplot as plt

data = np.loadtxt("z.dat")
data = data[np.newaxis, :]

# Determine array dimensions
N = data.shape[1]
Nhalf = int(np.fix(N/2))
is_odd = 0 if Nhalf*2==N else 1

# Do fft
freqs = np.concatenate([np.arange(0, Nhalf+1), np.arange(-Nhalf+1-is_odd, 0)]) / N
freqs = freqs[np.newaxis,:]
fdata = fftpack.fft(data)

# Compute all frequency domain Gaussians in one matrix
invfreqs = np.divide(1, freqs[:,1:Nhalf+1]).T
W = 2 * np.pi * np.multiply(matlib.repmat(freqs,Nhalf,1) , matlib.repmat(invfreqs,1,N))
G = np.exp((- W**2) / 2) # Gaussian in frequency domain

# Compute the Toeplitz matrix with the shifted frequency-domain data
HW = linalg.toeplitz(fdata[:, 0:Nhalf+1].T, fdata)
print(HW.shape)

# Exclude the first row, corresponding to zero frequency
HW = HW[1:Nhalf+1, :]

# Compute Stockwell Transform
ST = fftpack.ifft(np.multiply(HW, G), axis=1)

# Add the zero frequency row
ST0 = np.mean(data, axis=1) * np.ones(shape=(1, N))
ST = np.vstack((ST0, ST))

# Plot result
fig, ax1 = plt.subplots(1, 1, figsize=(10, 5))
im = ax1.pcolor(np.abs(ST))
ax1.set_xlabel("Time sample")
ax1.set_ylabel("Frequency")
ax1.set_title("S-Transform using my Python translation")
plt.colorbar(im)
plt.show()
plt.close()

问题根源 #

1. FFT/IFT归一化差异：Matlab的ifft默认会对结果做归一化（除以N），而SciPy的fftpack.ifft默认不归一化，这会直接导致幅度偏差，同时影响高斯窗的作用效果。
2. Toeplitz矩阵构造错误：Matlab的toeplitz(c,r)要求c为列向量、r为行向量，而原Python代码中对fdata的转置处理有误，导致矩阵结构与Matlab不一致。
3. 频率轴逻辑偏差：Matlab与Python的FFT负频率索引生成逻辑有细微差异，未对齐会导致高斯窗的计算错位。
4. 绘图轴方向差异：Matlab的axis xy会将y轴原点置于底部，而Pythonpcolor默认y轴原点在顶部，视觉上造成扭曲。

修正后的Python代码 #

import numpy as np
from scipy import fftpack, linalg
import matplotlib.pyplot as plt

# 读取数据并转为行向量
data = np.loadtxt("z.dat")
data = data.reshape(1, -1)

N = data.shape[1]
Nhalf = int(np.floor(N / 2))
is_odd = N % 2

# 生成与Matlab对齐的频率轴
if is_odd:
    freqs = np.concatenate([np.arange(0, Nhalf+1), np.arange(-Nhalf, 0)]) / N
else:
    freqs = np.concatenate([np.arange(0, Nhalf), np.arange(-Nhalf, 0)]) / N
freqs = freqs.reshape(1, -1)

# 执行FFT（与Matlab默认行为一致，不归一化）
fdata = fftpack.fft(data)

# 构造高斯窗矩阵，替换matlib.repmat避免依赖
inv_freqs = 1 / freqs[:, 1:Nhalf+1].T
freqs_rep = np.tile(freqs, (Nhalf, 1))
inv_freqs_rep = np.tile(inv_freqs, (1, N))
W = 2 * np.pi * freqs_rep * inv_freqs_rep
G = np.exp(-(W ** 2) / 2)

# 修正Toeplitz矩阵构造，匹配Matlab参数要求
c = fdata[:, 0:Nhalf+1].T
r = fdata[0, :]
HW = linalg.toeplitz(c, r)

# 移除0频率对应的行
HW = HW[1:Nhalf+1, :]

# 执行IFFT并添加归一化，对齐Matlab ifft行为
ST = fftpack.ifft(HW * G, axis=1, norm="ortho")

# 添加0频率行
ST0 = np.mean(data, axis=1).reshape(1, -1) * np.ones((1, N))
ST = np.vstack((ST0, ST))

# 绘图：对齐Matlab的axis xy显示逻辑
fig, ax1 = plt.subplots(1, 1, figsize=(10, 5))
im = ax1.imshow(np.abs(ST), aspect='auto', extent=[0, N, 0, ST.shape[0]])
ax1.set_xlabel("Time sample")
ax1.set_ylabel("Frequency")
ax1.set_title("S-Transform (Corrected Python Version)")
ax1.invert_yaxis()  # 将y轴原点移至底部，匹配Matlab显示
plt.colorbar(im)
plt.show()

关键修正说明 #

• 归一化对齐：在ifft中添加norm="ortho"，确保结果幅度与Matlab一致。
• Toeplitz矩阵修复：明确区分列向量参数c和行向量参数r，保证矩阵结构完全匹配Matlab。
• 频率轴校准：根据数据长度的奇偶性调整负频率生成逻辑，对齐Matlab的FFT频率索引。
• 绘图轴修正：使用invert_yaxis()将y轴方向反转，匹配Matlabaxis xy的视觉效果，消除扭曲。

内容的提问来源于stack exchange，提问作者user169511