speechproc.m

function speechproc()

    % 定义常数
    FL = 80;                % 帧长
    WL = 240;               % 窗长
    P = 10;                 % 预测系数个数
    s = readspeech('voice.pcm',100000);             % 载入语音s
    figure(1);
    subplot(2,1,1);
    plot(s);
    xlabel('origin');
    L = length(s);          % 读入语音长度
    FN = floor(L/FL)-2;     % 计算帧数
    % 预测和重建滤波器
    exc = zeros(L,1);       % 激励信号（预测误差）
    zi_pre = zeros(P,1);    % 预测滤波器的状态
    s_rec = zeros(L,1);     % 重建语音
    zi_rec = zeros(P,1);
    % 合成滤波器
    exc_syn = zeros(L,1);   % 合成的激励信号（脉冲串）
    s_syn = zeros(L,1);     % 合成语音
    % 变调不变速滤波器
    exc_syn_t = zeros(L,1);   % 合成的激励信号（脉冲串）
    s_syn_t = zeros(L,1);     % 合成语音
    % 变速不变调滤波器（假设速度减慢一倍）
    exc_syn_v = zeros(2*L,1);   % 合成的激励信号（脉冲串）
    s_syn_v = zeros(2*L,1);     % 合成语音
    FL_v = 2 * FL; 

    hw = hamming(WL);       % 汉明窗
    
    % 依次处理每帧语音
    for n = 3:FN

        % 计算预测系数（不需要掌握）
        s_w = s(n*FL-WL+1:n*FL).*hw;    %汉明窗加权后的语音
        [A, E] = lpc(s_w, P);            %用线性预测法计算P个预测系数
                                        % A是预测系数，E会被用来计算合成激励的能量

        A_rot = increase_150Hz(A);
        if n ==27
        % (3) 在此位置写程序，观察预测系统的零极点图
            % zplane(1, A);
            % figure(10);
            % freqz(1, A, 200);
            % figure(11);
            % freqz(1, A_rot, 200);
        end
        
        s_f = s((n-1)*FL+1:n*FL);       % 本帧语音，下面就要对它做处理

        % (4) 在此位置写程序，用filter函数s_f计算激励，注意保持滤波器状态
        [exc_tmp, zi_pre] = filter(A, 1, s_f, zi_pre);
        exc((n-1)*FL+1:n*FL) = exc_tmp; % 将你计算得到的激励写在这里
    
        % (5) 在此位置写程序，用filter函数和exc重建语音，注意保持滤波器状态
        [s_rec_tmp, zi_rec] = filter(1, A, exc_tmp, zi_rec);
        s_rec((n-1)*FL+1:n*FL) = s_rec_tmp; % 将你计算得到的重建语音写在这里
        
        % 注意下面只有在得到exc后才会计算正确
        s_Pitch = exc(n*FL-222:n*FL);
        PT = findpitch(s_Pitch);    % 计算基音周期PT（不要求掌握）
        G = sqrt(E*PT);           % 计算合成激励的能量G（不要求掌握）

        
        % (10) 在此位置写程序，生成合成激励，并用激励和filter函数产生合成语音
        % tmp = 0: (FL - 1);
        % exc_syn_tmp = double((mod(tmp, PT) == 0));
        % s_syn_tmp = filter(1, A, exc_syn_tmp);
        % s_syn_tmp = G * s_syn_tmp;

        % exc_syn((n-1)*FL+1:n*FL) = exc_syn_tmp; % 将你计算得到的合成激励写在这里
        % s_syn((n-1)*FL+1:n*FL) = s_syn_tmp; % 将你计算得到的合成语音写在这里

        % (11) 不改变基音周期和预测系数，将合成激励的长度增加一倍，再作为filter
        % 的输入得到新的合成语音，听一听是不是速度变慢了，但音调没有变。
        tmp = 0: (FL_v - 1);
        exc_syn_v_tmp = double((mod(tmp, PT) == 0));
        s_syn_v_tmp = filter(1, A, exc_syn_v_tmp);
        s_syn_v_tmp = G * s_syn_v_tmp;

        exc_syn_v((n-1)*FL_v+1:n*FL_v) = exc_syn_v_tmp; % 将你计算得到的加长合成激励写在这里
        s_syn_v((n-1)*FL_v+1:n*FL_v) = s_syn_v_tmp;  % 将你计算得到的加长合成语音写在这里

        
        % (13) 将基音周期减小一半，将共振峰频率增加150Hz，重新合成语音，听听是啥感受～
        % figure(10);
        % freqz(1, A, 200);
        % figure(11);
        % freqz(1, A_rot, 200);
        tmp = 0: (FL - 1);
        exc_syn_t_tmp = double( (mode(tmp, floor(PT/2)) == 0) );
        s_syn_t_tmp = filter(1, A_rot, exc_syn_t_tmp);
        s_syn_t_tmp = G * s_syn_t_tmp;

        exc_syn_t((n-1)*FL+1:n*FL) = exc_syn_t_tmp; % 将你计算得到的变调合成激励写在这里
        s_syn_t((n-1)*FL+1:n*FL) = s_syn_t_tmp; % 将你计算得到的变调合成语音写在这里
        
    end

    % (6) 在此位置写程序，听一听 s ，exc 和 s_rec 有何区别，解释这种区别
    sound(s);
    % sound(exc);
    % sound(s_rec);
    % sound(exc_syn);
    % figure(1);
    % stem(exc_syn);
    % sound(s_syn);
    subplot(2,1,2);
    plot(s_syn_t);
    xlabel('t-synthesis');
    sound(s_syn_t);
    

    % 保存所有文件
    writespeech('exc.pcm',exc);
    writespeech('rec.pcm',s_rec);
    writespeech('exc_syn.pcm',exc_syn);
    writespeech('syn.pcm',s_syn);
    writespeech('exc_syn_t.pcm',exc_syn_t);
    writespeech('syn_t.pcm',s_syn_t);
    writespeech('exc_syn_v.pcm',exc_syn_v);
    writespeech('syn_v.pcm',s_syn_v);
return

% 从PCM文件中读入语音
function s = readspeech(filename, L)
    fid = fopen(filename, 'r');
    s = fread(fid, L, 'int16');
    fclose(fid);
return

% 写语音到PCM文件中
function writespeech(filename,s)
    fid = fopen(filename,'w');
    fwrite(fid, s, 'int16');
    fclose(fid);
return

% 计算一段语音的基音周期，不要求掌握
function PT = findpitch(s)
    [B, A] = butter(5, 700/4000);
    s = filter(B,A,s);
    R = zeros(143,1);
    for k=1:143
        R(k) = s(144:223)'*s(144-k:223-k);
    end
    [R1,T1] = max(R(80:143));
    T1 = T1 + 79;
    R1 = R1/(norm(s(144-T1:223-T1))+1);
    [R2,T2] = max(R(40:79));
    T2 = T2 + 39;
    R2 = R2/(norm(s(144-T2:223-T2))+1);
    [R3,T3] = max(R(20:39));
    T3 = T3 + 19;
    R3 = R3/(norm(s(144-T3:223-T3))+1);
    Top = T1;
    Rop = R1;
    if R2 >= 0.85*Rop
        Rop = R2;
        Top = T2;
    end
    if R3 > 0.85*Rop
        Rop = R3;
        Top = T3;
    end
    PT = Top;
return