open_Call_b2g();open_Call_g2b();

% --- Executes on button press in open_pushbutton.
function open_Call_b2g()
% hObject    handle to open_pushbutton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
[filename pathname]=uigetfile({'*.wav','ALL FILES(*.*)'},'选择声音文件');
if isequal([filename pathname],[0,0])
    return;
end
strname=[pathname filename];%选择的声音文件路径和文件名
[data Fs]=wavread(strname);%data表示声音数据 Fs表示频率
handles.y=data;
handles.Fs=Fs;
handles.yoriginal=data;
handles.foriginal=Fs;
%guidata(hObject,handles);
b2g_Call(handles);

% --- Executes on button press in b2g_pushbutton.
function b2g_Call(handles)
% hObject    handle to b2g_pushbutton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
      %读取音频信息（双声道，16位，频率44100Hz）
      % 定义常数
    FL = 80;                % 帧长
    WL = 240;               % 窗长
    P = 10;                 % 预测系数个数
    data=handles.y;fs=handles.Fs;     % 载入语音数据
	data= data/max(data);	%归一化
    L = length(data);          % 读入语音长度
    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);     % 合成语音
	last_syn = 0;   %存储上一个（或多个）段的最后一个脉冲的下标
	zi_syn = zeros(P,1);   % 合成滤波器的状态
    
	% 变调不变速滤波器
    exc_syn_t = zeros(L,1);   % 合成的激励信号（脉冲串）
    s_syn_t = zeros(L,1);     % 合成语音
	last_syn_t = 0;   %存储上一个（或多个）段的最后一个脉冲的下标
	zi_syn_t = zeros(P,1);   % 合成滤波器的状态
    
	% 变速不变调滤波器（假设速度减慢一倍）
	
    hw = hamming(WL);       % 汉明窗
    
    % 依次处理每帧语音
    for n = 3:FN

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

      
        
        s_f = data((n-1)*FL+1:n*FL);       % 本帧语音，下面就要对它做处理

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

        % (5) 用filter函数和exc重建语音，注意保持滤波器状态
		[s_rec1,zi_rec] = filter(1,A,exc1,zi_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（不要求掌握）

        		
		%方法3：本段激励只能修改本段长度
		tempn_syn = [1:n*FL-last_syn]';
		exc_syn1 = zeros(length(tempn_syn),1);
		exc_syn1(mod(tempn_syn,PT)==0) = G; %某一段算出的脉冲
		exc_syn1 = exc_syn1((n-1)*FL-last_syn+1:n*FL-last_syn);
		[s_syn1,zi_syn] = filter(1,A,exc_syn1,zi_syn);
		exc_syn((n-1)*FL+1:n*FL) =  exc_syn1;   %计算得到的合成激励
		s_syn((n-1)*FL+1:n*FL) = s_syn1;   %计算得到的合成语音
		last_syn = last_syn+PT*floor((n*FL-last_syn)/PT);
 		
       
        
        % (13) 将基音周期减小一半，将共振峰频率增加150Hz，重新合成语音，听听是啥感受～
		PT1 =floor(PT/2);   %减小基音周期
        poles = roots(A);
		deltaOMG =150*2*pi/fs;
		for p=1:10   %增加共振峰频率，实轴上方的极点逆时针转，下方顺时针转
			if imag(poles(p))>0 poles(p) = poles(p)*exp(j*deltaOMG);
			elseif imag(poles(p))<0 poles(p) = poles(p)*exp(-j*deltaOMG);
			end
		end
		A1=poly(poles);
	
		
		tempn_syn_t = [1:n*FL-last_syn_t]';
		exc_syn1_t = zeros(length(tempn_syn_t),1);
		exc_syn1_t(mod(tempn_syn_t,PT1)==0) = G; %某一段算出的脉冲
		exc_syn1_t = exc_syn1_t((n-1)*FL-last_syn_t+1:n*FL-last_syn_t);
		[s_syn1_t,zi_syn_t] = filter(1,A1,exc_syn1_t,zi_syn_t);
		exc_syn_t((n-1)*FL+1:n*FL) =  exc_syn1_t;   %计算得到的合成激励
		s_syn_t((n-1)*FL+1:n*FL) = s_syn1_t;   %计算得到的合成语音
		last_syn_t = last_syn_t+PT1*floor((n*FL-last_syn_t)/PT1);
        
    end   %循环结束，获得变换后语音 s_syn_t
% 
%      plot(handles.axes3,s_syn_t),
%      set(handles.axes3,'Xgrid','on');
%      set(handles.axes3,'Ygrid','on');
%      xlabel(handles.axes3,'数据序列');
%      ylabel(handles.axes3,'频率');
%      title(handles.axes3,'变音后的时域图'),XLim([0,length(s_syn_t)]);	
% 
     subplot(1,3,3);plot(s_syn_t),
     grid on;
     xlabel('数据序列');
     ylabel('频率');
     title('变音后的时域图'),XLim([0,length(s_syn_t)]);	

     handles.y=s_syn_t;
%      guidata(hObject,handles);
     Fs=handles.foriginal;  
     N=length(handles.y);
    y1=fft(handles.y,N);
    f1=0:Fs/N:Fs*(N-1)/N;
%     plot(handles.axes1,f1,handles.yoriginal);
%     title(handles.axes1,'原始时域图');
%     xlabel(handles.axes1,'时间');
%     ylabel(handles.axes1,'频率');
%     set(handles.axes1,'Xgrid','on');
%     set(handles.axes1,'Ygrid','on');
%
    subplot(1,3,1);plot(f1,handles.yoriginal);
    title('原始时域图');
    xlabel('时间');
    ylabel('频率');
    grid on;

    
%     plot(handles.axes2,f1(1:N/2),y1(1:N/2));
%     set(handles.axes2,'Xgrid','on');
%     set(handles.axes2,'Ygrid','on');
%     title(handles.axes2,'频谱图');
%     xlabel( handles.axes2,'频率');
%     ylabel( handles.axes2,'幅度');
%
    subplot(1,3,2);plot(f1(1:N/2),y1(1:N/2));
    grid on;
    title('频谱图');
    xlabel('频率');
    ylabel('幅度');

   pause(3);
   sound(handles.yoriginal);
   sound(handles.y);

  % 计算一段语音的基音周期
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

   

% --- Executes on button press in open_pushbutton.
function open_Call_g2b()
% hObject    handle to open_pushbutton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
[filename pathname]=uigetfile({'*.wav','ALL FILES(*.*)'},'选择声音文件');
if isequal([filename pathname],[0,0])
    return;
end
strname=[pathname filename];%选择的声音文件路径和文件名
[data Fs]=wavread(strname);%data表示声音数据 Fs表示频率
handles.y=data;
handles.Fs=Fs;
handles.yoriginal=data;
handles.foriginal=Fs;
%guidata(hObject,handles);
g2b_Call(handles);

% --- Executes on button press in g2b_pushbutton.
function g2b_Call(handles)
% hObject    handle to g2b_pushbutton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
 y=handles.y;
 fs=handles.Fs;%读取音频信息（双声道，16位，频率44100Hz）
N=length(y)
% f=0:fs/N:fs*(N-1)/N;
f=[0:N-1]*fs/N;
Y=fft(y,N);    %进行傅立叶变换=>Y
%plot(handles.axes2,f(1:N/2),Y(1:N/2));
% title(handles.axes2,'声音信号的频谱');
% xlabel(handles.axes2,'频率');
% ylabel(handles.axes2,'振幅');
subplot(1,3,2);plot(f(1:N/2),Y(1:N/2));
title('声音信号的频谱');
xlabel('频率');
ylabel('振幅');
grid on;
% f1=0:(fs*0.7)/N:(fs*0.7)*(N-1)/N;
f1=[0:N-1]*(fs*0.7)/N;%降低频率为原频率f的0.7
% syms t;
t=[0,9];
R=y*exp(2*pi*300*t);%改变声音数据y=>R
P=fft(R,N);%y =>R , 傅立叶变换=>P
Z=ifft(P);
% z=real(Z);
handles.y=y; % y根本就没有变化啊?
% plot(handles.axes3,f1(1:N/2),Z(1:N/2));
% title(handles.axes3,'变声后的时域图');
% xlabel(handles.axes3,'时间序列');
% ylabel(handles.axes3,'频率')
% set(handles.axes3,'Xgrid','on');
% set(handles.axes3,'Ygrid','on');
subplot(1,3,3);plot(f1(1:N/2),Z(1:N/2));
title('变声后的频谱');
xlabel('频率');
ylabel('幅值')
grid on;

 Fs=handles.foriginal;  
    N=length(handles.y);
    y1=fft(handles.y,N);
%     f1=0:Fs/N:Fs*(N-1)/N;
    f1=[0:N-1]*Fs/N;
%     plot(handles.axes1,f1,handles.yoriginal);
%     title(handles.axes1,'原始时域图');
%     xlabel(handles.axes1,'时间');
%     ylabel(handles.axes1,'频率');
%     set(handles.axes1,'Xgrid','on');
%     set(handles.axes1,'Ygrid','on');
   subplot(1,3,1); plot(f1,handles.yoriginal);
    title('原始的频谱?');
    xlabel('频率');
    ylabel('幅值');
    grid on;
   
    
%     plot(handles.axes2,f1(1:N/2),y1(1:N/2));
%     set(handles.axes2,'Xgrid','on');
%     set(handles.axes2,'Ygrid','on');
%     title(handles.axes2,'频谱图');
%     xlabel( handles.axes2,'频率');
%     ylabel( handles.axes2,'幅度');
    figure;plot(f1(1:N/2),y1(1:N/2));
    grid on;
    title('频谱图');
    xlabel('频率');
    ylabel('幅度');

   pause(3);
%guidata(hObject,handles);

sound(handles.yoriginal);
sound(handles.y,handles.foriginal*0.7);