🎵 DSP 项目3 - 音频变调实验代码模板

📌 项目说明:
⚠️ 提交清单:
✅ 实验流程:
  1. 阅读实验手册《人声变调》(proj3_handout)
  2. 按顺序完成 ex0, ex2, ex3, ex5,填写标记为 ? 的代码
  3. 运行 test_Project3.p 检查代码
  4. 输入 9 生成加密成绩,将文件压缩为【学号_姓名.zip】提交
📄 ex0_all.m(主函数)需修改
function [y_interpolated, fs] = ex0_all(audio_path, alpha, Lwin, Ra, useMethod, beta)
% EX0_ALL - 根据指定参数对音频文件进行变调处理
%
% 输入参数:
%   audio_path  - 输入音频文件的路径
%   alpha       - 变调因子 (alpha > 1 升调, alpha < 1 降调)
%   Lwin        - 窗长
%   Ra          - 分析帧移
%   useMethod   - 1 使用方法一,2 使用方法二
%   beta        - 是否启用共振峰修正(可选,默认为1)

% 处理可选参数
if nargin < 6
    beta = 1;
end

% 清空变量和关闭图形窗口
clc;
clearvars -except audio_path alpha Lwin Ra useMethod beta;
close all;

% 读取输入音频文件
[x, fs] = audioread(audio_path);

% 确保信号为列向量
if size(x, 2) > 1
    x = x(:, 1);
end
x = x(:);

% 对输入音频进行 STFT 频谱分析
Freq = ex1_analysis(x, fs, Lwin, Ra);

% 根据选择的方法对频谱进行处理
if useMethod == 1
    % 方法1:Phase Vocoder,改变合成帧移实现变速
    Freq_out = ex2_process_method1(Freq, Lwin, Ra, alpha);
    % TODO: 方法1的合成帧移 = ?
    synthesize_param = ?;
elseif useMethod == 2
    % 方法2:频谱帧操作
    Freq_out = ex2_process_method2(Freq, Lwin, Ra, alpha);
    % TODO: 方法2的合成帧移 = ?
    synthesize_param = ?;
else
    error('useMethod 参数必须为 1 或 2');
end

% 调整共振峰
Freq_out_tuned = ex3_formant(Freq_out, alpha);

% 将处理后的频谱通过 ISTFT 合成回时域信号
y_out = ex4_synthesis(Freq_out_tuned, Lwin, fs, synthesize_param);

% 调整输出信号长度,使其与原始信号长度一致
y_interpolated = ex5_reformation(y_out, length(x));

% 播放处理后的音频
sound(y_interpolated, fs);
end
📄 ex2_process_method1.m(方法一变调)需修改
function [Freq_out] = ex2_process_method1(Freq, Lwin, Ra, alpha)
% EX2_PROCESS_METHOD1 频谱处理(方法一:Phase Vocoder)
%
% 方法一原理:通过修改合成帧移Rs实现时间伸缩,同时修正相位保证频率连续性

% 获取频谱矩阵的大小
[num_freqs, num_frames] = size(Freq);

% 计算频率索引和理论角频率
k = (0:num_freqs-1)';
Omega = 2 * pi * k / Lwin;

% 计算合成帧移
Rs = Ra * alpha;

% 初始化相位变量
prev_phase = zeros(num_freqs, 1);
phase_out = zeros(num_freqs, 1);

% TODO: 初始化输出频谱矩阵
Freq_out = ?;

% 从第2帧开始处理
for idx = ?
    % 提取当前帧的幅度和相位
    Amp = abs(Freq(:, idx));
    Pha = angle(Freq(:, idx));

    % TODO: 计算相位增量(相位展开)
    % 公式: mod(当前相位 - 上一帧相位 - 理论增量 + π, 2π) - π
    delta_phase = ?;

    % 更新上一帧相位
    prev_phase = Pha;

    % TODO: 计算频率偏移量
    delta_w = ?;
    
    % TODO: 累积输出相位(使用真实频率 × 合成帧移)
    phase_out = ?;

    % TODO: 生成输出频谱(幅度 × 复指数相位)
    Freq_out(:, idx) = ?;
end
end
📄 ex2_process_method2.m(方法二变调)需修改
function [Freq_out] = ex2_process_method2(Freq, Lwin, Ra, alpha)
% EX2_PROCESS_METHOD2 频谱处理(方法二:频谱帧操作)
%
% 方法二原理:通过对频谱帧进行抽取或复制改变帧数,重建时使用原帧移Ra

% 获取频谱矩阵的大小
[num_freqs, num_frames] = size(Freq);

% 计算频率索引和理论角频率
k = (0:num_freqs-1)';
Omega = 2 * pi * k / Lwin;

% 初始化累积相位
phase_out = zeros(num_freqs, 1);

% TODO: 计算新的时间索引
% 例如:alpha=2时,帧数增加;alpha=0.5时,帧数减少
t = ?;

% 计算输出帧数
num_frames_out = length(t);

% TODO: 初始化输出频谱矩阵
Freq_out = ?;

% 初始化帧计数器
itr = ?;

% 遍历每个新的时间索引
for idx = t
    % TODO: 计算当前帧索引
    col = ?;
    
    % 边界检查
    if col >= num_frames
        col = num_frames - 1;
    end
    if col < 1
        col = 1;
    end
    
    % 计算插值系数
    delta = idx - col;
    
    % 获取相邻帧幅值并进行线性插值
    p0Amp = abs(Freq(:, col));
    p1Amp = abs(Freq(:, min(col + 1, num_frames)));
    Amp = (1 - delta) * p0Amp + delta * p1Amp;
    
    % 获取相邻帧相位
    p0Pha = angle(Freq(:, col));
    p1Pha = angle(Freq(:, min(col + 1, num_frames)));
    
    % 计算相位增量(相位展开)
    delta_phase = mod(p1Pha - p0Pha - Omega * Ra + pi, 2 * pi) - pi;
    
    % 累积相位
    phase_out = phase_out + Omega * Ra + delta_phase;
    
    % TODO: 生成输出频谱
    Freq_out(:, itr) = ?;
    itr = itr + 1;
end
end
📄 ex3_formant.m(共振峰修正)需修改
function new_Freq = ex3_formant(Freq_out, alpha)
% EX3_FORMANT 对输入的频谱进行共振峰修正
%
% 共振峰修正原理:
% 1. 提取频谱包络(代表共振峰位置)
% 2. 计算残差 = 频谱 / 包络
% 3. 根据变调因子调整包络
% 4. 输出 = 残差 × 调整后的包络

[num_freqs, num_frames] = size(Freq_out);
new_Freq = zeros(size(Freq_out));
new_Freq(:, 1) = Freq_out(:, 1);

for idx = 2:num_frames
    % 提取幅度谱和相位谱
    Amp = abs(Freq_out(:, idx));
    phase_out = angle(Freq_out(:, idx));
    Amp_len = length(Amp);

    % 找到幅度谱的峰值作为共振峰估计
    [pks, loc] = findpeaks(Amp, 'MinPeakDistance', max(1, floor(Amp_len/20)));
    
    if isempty(pks)
        env = Amp;
    else
        % 对峰值进行插值得到包络
        loc_extended = [1; loc(:); Amp_len];
        pks_extended = [Amp(1); pks(:); Amp(end)];
        env = interp1(loc_extended, pks_extended, (1:Amp_len)', 'linear', 'extrap');
    end

    env(isnan(env)) = eps;
    env = max(env, eps);

    % 计算残差谱
    indices = env > eps;
    Amp_res = Amp;
    Amp_res(indices) = Amp(indices) ./ env(indices);

    % 根据alpha调整包络
    if alpha > 1
        cut = max(1, fix(Amp_len / alpha));
        spec_donor = zeros(size(Amp));
        old_indices = linspace(1, length(env), cut);
        new_indices = 1:cut;
        spec_donor(new_indices) = interp1(1:length(env), env, old_indices, 'linear', 'extrap');
    else
        cut = max(1, fix(length(env) * alpha));
        old_indices = 1:cut;
        new_indices = linspace(1, cut, Amp_len);
        spec_donor = interp1(old_indices, env(1:cut), new_indices, 'linear', 'extrap');
        spec_donor = spec_donor(:);
    end
    
    spec_donor(isnan(spec_donor)) = eps;
    spec_donor = max(spec_donor, 0);

    % 合成新的幅度谱
    new_Amp = Amp_res .* spec_donor;

    % 低通滤波
    filter_order = min(8, floor(Amp_len/4));
    cutoff_freq = min(0.99, 300/Amp_len);
    
    if cutoff_freq > 0 && cutoff_freq < 1 && filter_order > 0
        % TODO: 设计巴特沃斯低通滤波器
        [b, a] = ?;
        % TODO: 计算滤波器频率响应
        [h, ~] = ?;
        new_Amp = new_Amp .* abs(h);
    end

    % TODO: 生成新的复数频谱
    new_Freq(:, idx) = ?;
end
end
📄 ex5_reformation.m(插值重建)需修改
function [y_interpolated] = ex5_reformation(y, lx)
% EX5_REFORMATION 对信号进行插值重建
%
% 将变速后的信号恢复到原始长度,实现"变调不变速"

% 确保输入为实数列向量
y = real(y(:));
ly = length(y);

if ly == lx
    y_interpolated = y;
else
    % TODO: 使用interp1进行线性插值
    x_old = 1:ly;
    x_new = linspace(1, ly, lx);
    y_interpolated = ?;
end

y_interpolated = y_interpolated(:);

% TODO: 获取最小值和最大值
min_val = ?;
max_val = ?;

% TODO: 归一化到 [-1, 1] 范围
if max_val > min_val
    y_interpolated = ?;
else
    y_interpolated = zeros(lx, 1);
end
end
🧪 完成后测试:

完成所有代码后,可使用以下命令测试:

% 示例1:方法一升调1.5倍(无共振峰修正)
[y_out, fs] = ex0_all('./audio/orig1.wav', 1.5, 1024, 256, 1, 0);
% 预期:输出长度与原始相同,音调升高

% 示例2:方法一升调1.5倍(含共振峰修正)
[y_out, fs] = ex0_all('./audio/orig1.wav', 1.5, 1024, 256, 1, 1);
% 预期:音调升高但音色自然,不会像"小黄人"