如何绘制音频的频响曲线

什么是频响曲线

频响曲线（Frequency Response Curve）是表示系统或设备在不同频率下的输出响应特性的曲线。它通常用于音频设备、电路、机械系统等领域。通过观察频响曲线，我们可以了解系统在不同频率下的行为特性，例如哪些频率被放大或衰减。

1.整体计算流程

首先简述一下计算的流程，其实主要就是短时傅里叶变换的过程。之后再做详细的解释。
1.参数设置： 定义傅里叶变换的窗长和帧移，以及重叠长度。
2.生成窗函数
3.分帧处理并计算频谱： 对音频信号进行分帧处理，逐帧应用窗函数，并计算每帧的快速傅里叶变换（FFT），然后累加每帧的幅值谱。
4.求平均频谱： 对累加的频谱进行平均，得到整体的频响曲线。
5.绘制频率响应曲线： 将平均的频谱转换为分贝值，并在对数坐标系下绘制频响曲线。

2.为什么要分帧加窗

2.1 傅里叶变换的局限性

传统傅里叶变换的假设：傅里叶变换假定信号是 周期性 且 时不变 的，也就是说，信号的频率成分在整个时间段内保持不变。
• 问题：
○ 实际中的信号往往是 非平稳信号（频率特性会随时间变化）。
○ 直接对整个信号进行傅里叶变换无法捕捉频率随时间变化的信息。

2.2 短时傅里叶变换的目标

短时傅里叶变换（STFT）通过将信号分成多个短时段（帧），对每一帧进行傅里叶变换，从而捕捉信号在时间和频率上的局部特性。
• 分帧原因：假设在短时间段内，信号可以被近似为平稳信号。
• 加窗原因：
○ 分帧直接截断信号会产生频谱泄漏（即频率混叠）,频谱泄露会使得某些频率成分混杂到其他频率上，从而影响频谱分析的准确性。
○ 窗函数（如汉宁窗）可以平滑帧的边界，减少高频干扰。

2.3 窗函数的作用

窗函数会影响信号的频谱分辨率和泄漏：
• 频谱分辨率：窗的宽度决定了频率轴上的分辨率。较宽的窗函数提供更高的频率分辨率。
• 频谱泄漏：窗函数的边界光滑程度减少了频率混叠问题。

3. 常见的窗函数

1. 矩形窗 (Rectangular Window)：
   $$w(n)=10\le n\le N-1$$

2. 汉宁窗 (Hann Window)：
   $$w(n)=0.5\left(1-\cos \left(\frac{2\pi n}{N-1}\right)\right)0\le n\le N-1$$

3. 汉明窗 (Hamming Window)：
   $$w(n)=0.54-0.46\cos \left(\frac{2\pi n}{N-1}\right)0\le n\le N-1$$

4. 布莱克曼窗 (Blackman Window)：
   $$w(n)=0.42-0.5\cos \left(\frac{2\pi n}{N-1}\right)+0.08\cos \left(\frac{4\pi n}{N-1}\right)0\le n\le N-1$$

5. 凯泽窗 (Kaiser Window)：
   $$w(n)=\frac{I_0\left(\beta \sqrt{1-{\left(\frac{2n}{N-1}-1\right)}^2}\right)}{I_0(\beta )}0\le n\le N-1$$
   其中，( I_0 ) 是修正贝塞尔函数的零阶。

4. 处理代码(matlab)

% 读取音频文件
[audio_data, fs] = audioread('linsweep_96k.wav'); % 替换为你的音频文件路径，将音频数据读入 workspace

% 转换为单声道
if size(audio_data, 2) > 1
    audio_data = mean(audio_data, 2); % 如果是立体声文件，将左右声道取平均值，转换为单声道信号
end

% 参数设置
FFT_LENGTH = 4096; % FFT 的长度，影响频率分辨率，越长分辨率越高
HOP_SIZE = 2048;   % 帧移，定义每帧起点间的间隔（帧移越小帧间重叠越大）
OVERLAP = FFT_LENGTH - HOP_SIZE; % 重叠长度，用于计算帧之间信号的重叠部分

% 生成汉宁窗
hanning_window = 0.5 * (1 - cos(2 * pi * (0:FFT_LENGTH-1)' / (FFT_LENGTH-1)));
% 汉宁窗是一种平滑加权函数，用于减少 FFT 中频谱泄漏现象

% 初始化频谱累加数组
averaged_output = zeros(FFT_LENGTH / 2, 1); 
% 存储累加后的频谱幅值，大小为 FFT 长度的一半（仅存储对称频谱的前半部分）

% 计算信号长度和总帧数
signal_length = length(audio_data); % 获取音频信号的总采样点数
num_frames = floor((signal_length - OVERLAP) / HOP_SIZE); 
% 计算总帧数，信号的长度减去重叠部分后，再根据帧移分段

% 分帧处理并计算频谱
for frame = 1:num_frames
    % 当前帧的起始和结束索引
    start_idx = (frame - 1) * HOP_SIZE + 1; % 当前帧在信号中的起始位置
    end_idx = start_idx + FFT_LENGTH - 1;  % 当前帧在信号中的结束位置

    % 提取当前帧的信号
    if end_idx <= signal_length
        frame_signal = audio_data(start_idx:end_idx); 
        % 如果帧在信号范围内，直接取出这一帧的信号数据
    else
        frame_signal = zeros(FFT_LENGTH, 1); 
        % 如果帧超出信号范围，超出部分补 0
        valid_len = signal_length - start_idx + 1; % 计算有效信号长度
        frame_signal(1:valid_len) = audio_data(start_idx:signal_length);
        % 用信号的有效部分填充当前帧，剩余部分补 0
    end

    % 加窗
    frame_signal = frame_signal .* hanning_window; 
    % 对当前帧信号应用汉宁窗，以减少频谱泄漏

    % 计算 FFT
    fft_output = fft(frame_signal, FFT_LENGTH); 
    % 对当前帧信号进行 FFT，得到频域表示
    fft_magnitude = abs(fft_output(1:FFT_LENGTH/2)); 
    % 提取频谱的前半部分幅度谱（因为后半部分是对称的）

    % 累加频谱
    averaged_output = averaged_output + fft_magnitude / FFT_LENGTH;
    % 将当前帧的幅度谱归一化后累加到 `averaged_output`
end

% 求平均频谱
averaged_output = averaged_output / num_frames; 
% 将累加的频谱除以帧数，得到信号的平均频谱，降低短时变化对结果的影响

% 计算频率轴
f = (0:FFT_LENGTH/2-1) * fs / FFT_LENGTH; 
% 根据 FFT 的长度和采样率计算频率轴，用于绘图

% 绘制频率响应曲线
figure;
semilogx(f, 20*log10(averaged_output + eps)); 
% 使用对数频率轴绘制频谱曲线，并将幅值转为分贝（加 `eps` 避免对 0 取 log）
xlabel('Frequency (Hz)'); % 标注横轴为频率
ylabel('Magnitude (dB)'); % 标注纵轴为幅值（分贝）
title('Frequency Response'); % 设置图表标题
grid on; % 打开网格线，便于观察

% 显示频率范围限制（例如 20 Hz 到 20 kHz）
xlim([20 20000]); % 限制横轴范围到人耳可听频率范围

4.1为什么需要累加频谱？

% 累加频谱
averaged_output = averaged_output + fft_magnitude / FFT_LENGTH;

对每一帧的 FFT 幅度谱进行归一化（通过 fft_magnitude / FFT_LENGTH），然后将归一化结果逐帧累加到 averaged_output 中。
音频信号通常是非平稳的，其频率特性可能在不同时间段有所变化。通过逐帧累加，能够整合多个时间窗口的频谱信息，从而反映信号的整体频率特性。
如果不累加，最后就只有一帧的信息，如下图：

之后再除以累加的帧长，得到平均响度

averaged_output = fft_magnitude / num_frames;