Python音频处理-频谱图实现

from IPython.display import Audio
from torchaudio.utils import download_asset
from librosa import power_to_db
import matplotlib.pyplot as plt
import numpy as np
import torchaudio

# 基本参数
n_fft = 512
hop_length = 512

# 示例音频
sample_speech = download_asset(
    "tutorial-assets/Lab41-SRI-VOiCES-src-sp0307-ch127535-sg0042.wav"
)

# 提取频谱
waveform, sample_rate = torchaudio.load(sample_speech)
print(waveform.shape)
print(sample_rate)

spectrogram = torchaudio.transforms.Spectrogram(
    n_fft=n_fft,
    hop_length=hop_length
)
spec = spectrogram(waveform)
spec.shape

# 绘制波形和频谱图
_, axs = plt.subplots(2, 1)
axs[0].plot(waveform[0])
axs[1].imshow(
    power_to_db(spec[0]), 
    origin='lower', 
    aspect='auto', 
    interpolation='nearest'
)
plt.tight_layout()
plt.show()

# 手动实现频谱计算
def my_spectrogram(signal, n_fft=512, hop_length=512):
    # 数组转换
    signal = np.asarray(signal)
    # 计算帧数
    signal_length = len(signal)
    num_frames = int(np.ceil(signal_length / hop_length))
    # 扩展信号
    pad_length = (num_frames - 1) * hop_length + n_fft - signal_length
    signal = np.pad(
        signal,
        (0, pad_length),
        mode='constant'
    )
    # 分帧
    frames = np.lib.stride_tricks.sliding_window_view(
        signal,
        window_shape=n_fft
    )[::hop_length]
    # 加窗
    window = np.hanning(n_fft)
    frames = frames * window
    # 快速傅里叶变换
    fft_result = np.fft.fft(frames, axis=1)
    magnitude_spectrum = np.abs(fft_result) ** 2
    # 对称处理
    magnitude_spectrum = magnitude_spectrum[:, : n_fft//2 + 1]
    return magnitude_spectrum

# 计算频谱
my_spec = my_spectrogram(waveform[0].numpy())
my_spec.shape

# 绘制波形和频谱图
_, axs = plt.subplots(2, 1)
axs[0].imshow(
    power_to_db(spec[0]), 
    origin='lower', 
    aspect='auto', 
    interpolation='nearest'
)
axs[1].imshow(
    power_to_db(my_spec.T), 
    origin='lower', 
    aspect='auto', 
    interpolation='nearest'
)
plt.tight_layout()
plt.show()

# 计算两个频谱差距
import torch
mse = torch.square(spec[0] - my_spec.T).mean().item()
mse