使用pytorch进行FFT和STFT

首先,我们定义一个波形,幅值分别为20和38,频率为2和13:
y = 20 sin ⁡ ( 2 π × 2 x ) + 38 sin ⁡ ( 2 π × 13 x ) y=20 \sin (2\pi \times 2x)+38\sin (2\pi \times 13x) y=20sin(2π×2x)+38sin(2π×13x)
采样频率为200Hz,采样时间为1s。由于Pytorch中没有类似于Numpy中numpy.pi的用法,所以我们先用Numpy计算函数,然后再存入Tensor。

import torch
import numpy as np
import matplotlib.pyplot as plt

Fs = 200 # Sampling frequency
T = 1 / Fs # period
t = 1 # time
L = Fs * t # signle length

x = np.arange(L)*T
y = 20*np.sin(2*np.pi*2*x) + 38*np.sin(2*np.pi*13*x)

plt.plot(x,y)
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.show()

波形图如下:
在这里插入图片描述
傅里叶变换是一种将信号从时域转换到频域的数学工具。它将一个连续或离散的信号分解成一系列不同频率的正弦和余弦波的叠加。使用pytorch对该函数进行快速傅里叶变换:

p = torch.fft.fft(torch.from_numpy(y))
p = 2 * torch.abs(p) / len(p)
p = p[:int(len(p) / 2)]
p[0] = 0.
freq_bin = np.arange(int(L/2))*Fs/L

plt.plot(freq_bin,p)
plt.xlabel('Frequency (Hz)')
plt.ylabel('|P1|')
plt.show()

频谱图如下,可以看到幅值和频率与函数y相同:
在这里插入图片描述
下面,我们再对函数进行短时傅里叶变换。短时傅里叶变换(Short-Time Fourier Transform,STFT)是一种信号处理技术,用于将非平稳信号分解成时间和频率的分量。它是傅里叶变换(Fourier Transform)的一种变体,在时间上加一定的窗函数,然后对每个窗函数内的信号进行傅里叶变换,从而得到时频域上的信号特征。STFT常用于音频、图像等信号处理领域,以提取时域和频域上的信息,以及信号的时频分布特性:

spectrum = torch.stft(input=torch.from_numpy(y), n_fft=200, return_complex=True)

此外,还可以用torchaudio.functional.spectrogram,也是返回STFT之后的频谱。

如果想直接画出频谱图可以使用matplotlib的plt.specgram方法:

plt.specgram(y, Fs=Fs)
plt.xlabel('Time (s)')
plt.ylabel('Frequency (Hz)')
plt.show()

效果如下:
在这里插入图片描述

  • 1
    点赞
  • 8
    收藏
    觉得还不错? 一键收藏
  • 1
    评论
您好!我了解您的问题,以下是使用 PyTorch 进行音频 STFT 的代码示例: 首先,您需要导入所需库: ``` python import torch from torch.nn.functional import pad # 如果您希望可视化输出,请安装 matplotlib 库 import matplotlib.pyplot as plt ``` 接下来,您需要定义一个函数来计算短时傅里叶变换STFT): ``` python def stft(x, fft_size, hop_size, win_size=None, window='hann'): # 将输入转换为可用于torch.fft的复数张量 x = x.unsqueeze(-1) x_real = torch.cat([x, torch.zeros_like(x)], dim=-1) # 创建窗口函数(默认为汉宁窗) if win_size is None: win_size = fft_size win = torch.hann_window(win_size, device=x.device) if window == 'hann' else torch.ones(win_size, device=x.device) # 零填充输入以获得所需帧长度 pad_size = (fft_size - hop_size) // 2 x_real = pad(x_real, [pad_size, pad_size], mode='reflect') # 计算每个帧的STFT stft_frames = [] for i in range((x_real.shape[-2] - fft_size) // hop_size + 1): start = i * hop_size end = start + fft_size frame = x_real[..., start:end] # 加窗 frame *= win # 计算FFT和取模 frame = torch.fft.fft(frame, fft_size, dim=-2) frame = frame.real**2 + frame.imag**2 stft_frames.append(frame) # 将所有帧堆叠成一张张量 stft_tensor = torch.stack(stft_frames, dim=-1) return stft_tensor ``` 这个函数的输入是一个一维张量 x,表示原始音频信号。fft_size 是要计算的FFT大小,hop_size 是每个帧之间的跳跃量。如果未指定 win_size,则使用 fft_size,窗口函数可以是 'hann' 或 'rect',默认为 'hann'。输出是一张三维张量,其形状为(freq_bins,frames,1),其中 freq_bins 是 FFT 的时间长度,frames 是帧的数量。 例如,如果您有一个名为 audio 的一维张量,采样率为 16kHz,想要在窗口大小为 1024,每隔 256 个样本(即每 16ms)计算一个帧的情况下计算 256 点的FFT,您可以这样调用函数: ``` python sample_rate = 16000 window_size = 1024 hop_size = 256 freq_bins = window_size // 2 + 1 audio = torch.randn(sample_rate * 10) # 10秒的随机信号 stft_matrix = stft(audio, window_size, hop_size) # 可视化stft(仅供参考) plt.imshow(stft_matrix[:, :, 0].log2().flip(0), cmap='jet', extent=[0, 10, 0, freq_bins], aspect='auto') plt.xlabel('Time (s)') plt.ylabel('Frequency (Hz)') plt.show() ``` 这将计算 stft_matrix,其中 freq_bins 是 513,帧数为 (len(audio) - window_size) / hop_size + 1。 请注意,该函数实现了“重叠相加”的 STFT,这意味着每个帧之间存在一些重叠,以使输出保持平滑,并且第一个样本之前和最后一个样本之后采用“反射”零填充以避免边缘效应。 希望这可以回答您的问题!

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值