DCT 和 FFT 的对比：短时傅里叶变换与短时离散余弦变换的实现

在信号处理领域，离散余弦变换（DCT）和快速傅里叶变换（FFT）是两种重要的变换方法，广泛应用于频域分析、特征提取和数据压缩。本文将对 DCT 和 FFT 进行对比，并介绍短时傅里叶变换（STFT）和短时离散余弦变换（STDCT）的物理意义、参数意义、计算公式以及数据输入输出维度。

1. STFT 和 STDCT 的对比

1.1 短时傅里叶变换（STFT）

物理意义

STFT 对信号进行短时分析，将其分为短时间帧并对每帧应用傅里叶变换，以捕捉信号在时间和频率上的变化。

参数意义

• nperseg：每个分段的长度，决定频率分辨率。
• noverlap：每个分段之间的重叠样本数，影响时间分辨率。

计算公式

对于信号 $x[n]$，STFT 的计算公式为：

$$X[k,m]=\sum _{n=0}^{N-1}x[n]w[n-m]e^{-j\frac{2\pi }{N}kn}$$

其中 $w[n]$ 是窗函数。

数据输入输出维度

• 输入：一维信号，形状为 $N$，其中 $N$ 是信号的长度。
• 输出：二维频谱，形状为 $(N_f,N_t)$，其中 $N_f$ 是频率分量的数量，$N_t$ 是时间帧的数量。

1.2 短时离散余弦变换（STDCT）

物理意义

STDCT 对信号进行短时分析，将其分为短时间帧并对每帧应用离散余弦变换，有效表示信号的频谱特性，尤其是在信号能量集中在低频部分时。

参数意义

• frame_length：每个帧的长度，影响频率分辨率。
• hop_length：帧移，影响时间分辨率。

计算公式

对于信号 $x[n]$，STDCT 的计算公式为：

$$X[k,m]=\sum _{n=0}^{N-1}x[n]w[n-m]\cos \left(\frac{\pi }{N}\left(n+\frac{1}{2}\right)k\right)$$

数据输入输出维度

• 输入：一维信号，形状为 $N$，其中 $N$ 是信号的长度。
• 输出：二维频谱，形状为 $(N_d,N_t)$，其中 $N_d$ 是 DCT 系数的数量（通常为帧长度），$N_t$ 是时间帧的数量。

2. Python 代码实现

以下是一个完整的 Python 示例代码，展示如何生成随机信号，计算 STFT 和 STDCT，并重建信号。

import numpy as np
import matplotlib.pyplot as plt
import librosa
import librosa.display
from scipy.fftpack import dct, idct

# 1. 生成随机信号
fs = 16000 # 采样率
t = np.linspace(0, 1, fs, endpoint=False)  # 1秒钟的时间
signal = np.random.randn(fs)  # 生成标准正态分布的随机信号

# 2. 使用 librosa 计算 STFT
n_fft = 256  # FFT 的长度
hop_length = 128  # 帧移，默认情况下可以设置为 n_fft // 2
Zxx_librosa = librosa.stft(signal, n_fft=n_fft, hop_length=hop_length)
print(Zxx_librosa.shape)

# 计算频率和时间轴
frequencies_librosa = np.fft.rfftfreq(n_fft, d=1 / fs)  # 计算频率
times_librosa = np.arange(Zxx_librosa.shape[1]) * hop_length / fs  # 计算时间

print("Number of frames (times_librosa.shape):", times_librosa.shape[0])


# 计算 STDCT
def short_time_dct(signal, frame_length=256, hop_length=128):
    # 窗函数
    window = np.hanning(frame_length)
    # 计算帧数
    num_frames = (len(signal) - frame_length) // hop_length + 1
    # 初始化 DCT 结果
    dct_result = np.zeros((num_frames, frame_length))

    for i in range(num_frames):
        start = i * hop_length
        frame = signal[start:start + frame_length] * window
        dct_result[i, :] = dct(frame, norm='ortho')

    return dct_result


stdct_result = short_time_dct(signal)
print("STDCT shape:", stdct_result.shape)


# 3. 重建信号
def reconstruct_signal_from_stft(Zxx, hop_length):
    # 使用 librosa 的 istft
    reconstructed_signal = librosa.istft(Zxx, hop_length=hop_length)
    return reconstructed_signal


def reconstruct_signal_from_stdct(dct_result, hop_length, frame_length):
    num_frames = dct_result.shape[0]
    reconstructed_signal = np.zeros(num_frames * hop_length + frame_length)
    window = np.hanning(frame_length)

    for i in range(num_frames):
        start = i * hop_length
        frame = idct(dct_result[i, :], norm='ortho')
        reconstructed_signal[start:start + frame_length] += frame * window

    return reconstructed_signal


reconstructed_signal_stft = reconstruct_signal_from_stft(Zxx_librosa, hop_length=hop_length)
reconstructed_signal_stdct = reconstruct_signal_from_stdct(stdct_result, hop_length=128, frame_length=256)

# 4. 可视化
plt.figure(figsize=(12, 10))

# 原始信号
plt.subplot(4, 1, 1)
plt.plot(t, signal)
plt.title('Original Signal')
plt.xlabel('Time [s]')
plt.ylabel('Amplitude')

# STFT 结果
plt.subplot(4, 1, 2)
plt.pcolormesh(times_librosa, frequencies_librosa, np.abs(Zxx_librosa), shading='gouraud')
plt.title('STFT Magnitude using librosa')
plt.xlabel('Time [s]')
plt.ylabel('Frequency [Hz]')

# STDCT 结果
plt.subplot(4, 1, 3)
plt.imshow(np.abs(stdct_result.T), aspect='auto', origin='lower', extent=[0, t[-1], 0, 256])
plt.title('STDCT Magnitude')
plt.xlabel('Time [s]')
plt.ylabel('DCT Coefficients')

# 重建信号
plt.subplot(4, 1, 4)
plt.plot(t, reconstructed_signal_stft[:len(t)], label='Reconstructed from STFT', alpha=0.7)
plt.plot(t, reconstructed_signal_stdct[:len(t)], label='Reconstructed from STDCT', alpha=0.7)
plt.title('Reconstructed Signals')
plt.xlabel('Time [s]')
plt.ylabel('Amplitude')
plt.legend()

plt.tight_layout()
plt.show()