教学用Matlab语音处理程序

function speechproc（）
clear allclose all clc;
    % 定义常数
    FL = 80;                % 帧长
    WL = 240;               % 窗长
    P = 10;                 % 预测系数个数
    s = readspeech（‘voice.pcm‘100000）;             % 载入语音s
	s = s/max（s）;	%归一化
    L = length（s）;          % 读入语音长度
    FN = floor（L/FL）-2;     % 计算帧数
    
	% 预测和重建滤波器
    exc = zeros（L1）;       % 激励信号（预测误差）
    zi_pre = zeros（P1）;    % 预测滤波器的状态
    s_rec = zeros（L1）;     % 重建语音
    zi_rec = zeros（P1）;
    
	% 合成滤波器
    exc_syn = zeros（L1）;   % 合成的激励信号（脉冲串）
    s_syn = zeros（L1）;     % 合成语音
	last_syn = 0;   %存储上一个（或多个）段的最后一个脉冲的下标
	zi_syn = zeros（P1）;   % 合成滤波器的状态
    
	% 变调不变速滤波器
    exc_syn_t = zeros（L1）;   % 合成的激励信号（脉冲串）
    s_syn_t = zeros（L1）;     % 合成语音
	last_syn_t = 0;   %存储上一个（或多个）段的最后一个脉冲的下标
	zi_syn_t = zeros（P1）;   % 合成滤波器的状态
    
	% 变速不变调滤波器（假设速度减慢一倍）
	v=.5;
    exc_syn_v = zeros（v\L1）;   % 合成的激励信号（脉冲串）
    s_syn_v = zeros（v\L1）;     % 合成语音
	last_syn_v = 0;   %存储上一个（或多个）段的最后一个脉冲的下标
	zi_syn_v = zeros（P1）;   % 合成滤波器的状态
	
    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会被用来计算合成激励的能量

        if n == 27
        % （3） 观察预测系统的零极点图
            zplane（1A）;
        end
        
        s_f = s（（n-1）*FL+1:n*FL）;       % 本帧语音，下面就要对它做处理

        % （4） 用filter函数s_f计算激励，注意保持滤波器状态
		[exc1zi_pre] = filter（A1s_fzi_pre）;
        
        exc（（n-1）*FL+1:n*FL） = exc1; %计算得到的激励

        % （5） 用filter函数和exc重建语音，注意保持滤波器状态
		[s_rec1zi_rec] = filter（1Aexc1zi_rec）;
        
        s_rec（（n-1）*FL+1:n*FL） = s_rec1; %计算得到的重建语音

        % 注意下面只有在得到exc后才会计算正确
        s_Pitch = exc（n*FL-222:n*FL）;
        PT = findpitch（s_Pitch）;    % 计算基音周期PT（不要求掌握）
        G = sqrt（E*PT）;           % 计算合成激励的能量G（不要求掌握）

        
        % （10） 生成合成激励，并用激励和filter函数产生合成语音
		% 方法1：本段激励可以向右超出段长一个基音周期
% 		tempn_syn = [1:n*FL-last_syn+PT-1]‘;
% 		exc_syn1 = zeros（length（tempn_syn）1）;
% 		exc_syn1（mod（tempn_synPT）==0） = G; %某一段算出的脉冲
% 		[s_syn1zi_syn] = filter（1Aexc_syn1zi_syn）;
% 		exc_syn（last_syn+1:last_syn+length（tempn_syn）） = exc_syn1;
% 		s_syn（last_syn+1:last_syn+length（tempn_syn）） = s_syn1;	
% 		last_syn = last_syn+PT*floor（length（tempn_syn）/PT）;
       
		%方法2：本段激励可以向左侵占前面一段
% 		tempn_syn = [1:n*FL-last_syn]‘;
% 		exc_syn1 = zeros（length（tempn_syn）1）;
% 		exc_syn1（mod（tempn_synPT）==0） = G; %某一段算出的脉冲
% 		[s_syn1zi_syn] = filter（1Aexc_syn1zi_syn）;		
% 		exc_syn（last_syn+1:n*FL） = exc_syn1;
% 		s_syn（last_syn+1:n*FL） = s_syn1;	
% 		last_syn = last_syn+PT*floor（（n*FL-last_syn）/PT）;
		
		%方法3：本段激励只能修改本段长度
		tempn_syn = [1:n*FL-last_syn]‘;
		exc_syn1 = zeros（length（tempn_syn）1）;
		exc_syn1（mod（tempn_synPT）==0） = G; %某一段算出的脉冲
		exc_syn1 = exc_syn1（（n-1）*FL-last_syn+1:n*FL-last_syn）;
		[s_syn1zi_syn] = filter（1Aexc_syn1zi_syn）;
		exc_syn（（n-1）*FL+1:n*FL） =  exc_syn1;   %计算得到的合成激励
		s_syn（