%=================================================================
%             PSOLA算法實現程序
%=================================================================
%
%具體的算法在畢業論文《基于PSOLA的漢語文語轉換技術研究》P44中有詳細的介紹。
%
%需要注意的問題：目前對聲母（清音）還沒有修改，只是修改韻母（濁音），在必要的時候可以適當加以處理聲母。
%對聲母的處理只能調整音長和音強，也就是復制一段接在聲母后面以增長音長，乘以某個系數以提高音強。
%=================================================================
%2006年7月    郭鋒    版本 1.0
%=================================================================
clear all
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%產生4個聲調所對應的聲調函數，數據通過對聲調曲線的圖像來采樣獲得
                                                              %實驗證實：11個采樣點已經可以得到實驗要求。
t=0:0.1:1.0;
f1=[0.453 0.470 0.472 0.469 0.466 0.465 0.468 0.473 0.476 0.470 0.447];
f2=[0.011 0.081 0.10 0.128 0.163 0.236 0.305 0.384 0.460 0.520 0.544];
f3=[-0.185 -0.106 -0.151 -0.205 -0.2309 -0.20 -0.105 0.001 0.124 0.220 0.324];
f4=[0.463 0.531 0.504 0.438 0.363 0.267 0.184 0.120 0.046 0.001 -0.064];
f1=f1+0.6;  f2=f2+0.6;  f3=f3+0.6;  f4=f4+0.6; %往上提升0.6，把數據全部化為正數
kp1=polyfit(t,f1,4);         %擬和得到四個聲調的聲調函數，它的圖像形式就是調型曲線，分別用kp1～kp4來標示。
kp2=polyfit(t,f2,4);
kp3=polyfit(t,f3,4);
kp4=polyfit(t,f4,4);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%產生基音變化序列
[pitch,Fs]=wavread('a.wav');                                      %讀入一個基音周期，例如a的一個基音。
OriginPitchLength=length(pitch);
PithNum=40;                   %基音數目從此輸入，對應著音長。
z=linspace(0,1,PithNum);      
Fn=polyval(kp4,z);            %通過調型函數來獲得需要的聲調序列。
Fn=1./Fn;                     %計算基頻數列
KK=log(10+10.*Fn);            
scale=OriginPitchLength/KK(PithNum/2);    %scale表示單位數據代表的基音點數
ResultPitchLength=round(scale.*KK);       %獲得該聲調對應的基音周期序列，準備用PSOLA調整基音周期。
%以下四行注釋為早期的一些實驗結果，這里不用理會。
%chazhiSerial=round(PitchSerial-PitchSerial(15))    %獲得插值點數序列
%測試二聲的序列  最大基音和最小基音相差50個  169.6154－－122.5000
%三聲   最大基音和最小基音相差40個   183.7500－－147
%四聲   最大基音和最小基音相差60個   122.5000－－183.7500
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%開始基音同步疊加,改變聲調
a=pitch';
resultpitch=a;  %開始調整基音resultpitch。
reginpitch=a;   %reginpitch留作最后來和調整過的resultpitch進行比較。
K=zeros(1,PithNum);
for ii=1:PithNum
    if (OriginPitchLength==ResultPitchLength(ii))        %如果需要的基音數值NeedNum與原始基音數值OriginNum相等
        c=a;                                             %那么就不作修改，直接把原始基音a給存儲變量c。
    else                                                 %否則的話，就對原始基音進行先插值再抽取。
         K(ii)=lcm(OriginPitchLength,ResultPitchLength(ii));  %K為原始基音數值，例如5，和需要的基音數值，例如7，的最小公倍數。
         D1(ii)=K(ii)/OriginPitchLength;
         D2(ii)=K(ii)/ResultPitchLength(ii);
         b=zeros(1,K(ii));
          num=1;
         %%%%%%%%%%%%%%%%%%%%%%%%%%%                     插值
         for i=1:K(ii)
           if(rem(i,D1(ii))==1)          %D1(ii)整倍，就直接賦值，不作修改
             b(i)=a(num);
             num=num+1;
           elseif (num<=OriginPitchLength) %否則，進行線性插值，擴充點數。
             b(i)=b(i-1)+(a(num)-a(num-1))/(D1(ii)-1);
           elseif (num>OriginPitchLength)
            b(i)=a(num-1);       %如果是最后一個數據點，就直接幅值，不需要再插值。
           end
         end
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  按照需要的間隔進行數據抽取
         count=1;
         c=zeros(1,ResultPitchLength(ii));
        for j=1:K(ii)
         if (rem(j,D2(ii))==0)
           c(count)=b(j);
           count=count+1;
         end
        end
    end
     reginpitch=[reginpitch a];
     resultpitch=[resultpitch c];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%比較結果
[am,an]=size(reginpitch);
[cm,cn]=size(resultpitch);
%reginpitch=[a zeros(1,abs(cn-an))];
[qingyin,Fs]=wavread('t.wav');      %讀入清音音節
%%%%%%%%%%%%%%%%%%%%%%%%%用曲線進行衰減語音波形，使得整個濁音語音波形中間音強高，首尾低，方才符合人類自然語音。
tjidian=1/100;
xjidian=0.8/65;
t0=[3 12 21 30 56 76 85 93 97];
x0=[23 32 43 41 37 33 28 22 19];
T0=tjidian.*t0;
X0=xjidian.*x0+0.2;
s1=polyfit(T0,X0,5);      %采用的衰減曲線s1為開口向下的拋物線。
%z=linspace(0,1,100);
%Fn=polyval(p1,z);
%plot(z,Fn)          %可以觀測衰減曲線
k1=linspace(0,1,an);
k2=linspace(0,1,cn);
sjn1=polyval(s1,k1);   %獲得原始基音對應的衰減序列。
sjn2=polyval(s1,k2);   %獲得合成基音對應的衰減序列。
sum1=reginpitch.*sjn1;  %衰減原始基音序列
sum2=resultpitch.*sjn2; %衰減合成基音序列
outSpeechResult=[qingyin' sum2];
outSpeechOrign=[qingyin'  sum1];    %給outSpeechOrign乘以系數K就可以修改音強。
wavplay(outSpeechResult,Fs)         %讀出合成語音
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%把結果以圖形的形式顯示出來。
subplot(2,1,1)
plot(outSpeechOrign)
grid on
title('漢語“他（ta1）”的時域波形圖')
wavplay(outSpeechOrign,Fs)
subplot(2,1,2)
plot(outSpeechResult)
title('TD-PSOLA處理后變為“踏（ta4）”的時域波形圖')
grid on