今天看到群里有人討論這個問題，記錄一下。
主要內容轉自：http://www.cnblogs.com/welen/p/3782896.html
變調和變速原理
自然語音的產生可以簡化為圖2-1模型，激勵源出來的聲門波信號與聲道模型進行卷積，最后通過嘴唇輻射模型產生語音。其中，激勵源決定說話人的基頻的大小，即音調的高低。聲道模型反映“潤色”的頻譜信息，具體的講，共振峰決定了語義信息，諧波分布決定了音色，單位時間的音節數決定了語速。
圖2-1 語音產生模型
下面將根據語音產生模型來闡述變速變調的基本原理。
變速變調的改變可以包括變速不變調和變調不變速兩個部分。
語音變速不變調是指保持音調和語義保持不變，語速變快或變慢[28]。該過程表現為語譜圖在時間軸上如手風琴般壓縮或者擴展。那也就是說，基頻值幾乎不變，對應於音調不變；整個時間過程被壓縮或者擴展，聲門周期的數目減小或者增加，即聲道運動速率發生改變，語速也隨之變化。對應於語音產生模型，激勵和系統經歷與原始發音情況幾乎相同的狀態，但持續時間相比原來或長或短[29]。
嚴格地講，基頻和音調是兩個不同的概念，基頻是指聲帶振動的頻率，音調是指人類對基頻的主觀感知，但是兩者變化基本一致，即基頻越高，音調越高，基頻越低，音調越低，音調是由基頻決定的[30]。因此，語音變調不變速就是指改變說話人基頻的大小[44]，同時保持語速和語義不變，即保持短時頻譜包絡(共振峰的位置和帶寬)和時間過程基本不變[31]。對應於語音產生模型，變調改變了激勵源；聲道模型的共振峰參數幾乎不變，保證了語義和語速不變。
綜上所述，變速改變聲道運動速率，力求保持激勵源不變；變調改變激勵源，力求保持聲道的共振峰信息不變。但是聲源和聲道不是相互獨立的，在改變聲源時，必然也會非線性的影響聲道，同樣地，改變聲道時也會或多或少的影響聲源，兩者之間相互影響，相互作用。
變調不變速方法
變調的方法也可以分為三類：時域法、頻域法、參量法。
時域法中，Crochiere等人於1983年提出了重采樣的方法[42]，該方法是實現變速變調最簡單、最常用的方法之一。
假設重采樣因子為P/Q，其中，P為上采樣因子，Q為下采樣因子。上采樣過程就是往原始信號相鄰兩點間內插P-1個采樣點，這樣使得基音周期變為原來的P倍，頻譜壓縮為原來的1/P倍，時長變為原來的P倍，即基頻變為原來的1/P倍，音調降為原來的1/P倍，語速變為原來的1/P倍。
同樣地，下采樣過程就是每隔Q-1個點進行抽取，這樣會使得基音周期長度為原來的1/Q倍，頻譜擴展為原來的Q倍，時長變為原來的1/Q倍，即基頻變為原來的Q倍，音調升為原來的Q倍，語速變為原來的Q倍。
綜合上述兩個過程，通過P/Q倍的重采樣后，保持播放速率不變，重采樣語音語速和音調都變為原來的Q/P倍[43]。
為了實現變調不變速，可以通過各種變速不變調處理與重采樣相結合的方法[44]。如圖2-4所示，變速不變調處理使語速變為原來的P/Q倍，得到輸出信號y(n)，然后對y(n)進行P/Q倍重采樣處理，這樣就得到語速正常，音調變為原來Q/P倍的最終輸出語音z(n)。
頻域法中比較簡單的處理就是直接對信號頻譜進行插值或者抽取，實現各頻率分量的擴展或者壓縮。國內的研究者李力利、張曉蕊等人分別對頻域的插值和抽取的方法進行了研究和擴展，這種方法的缺點在於：內插會引入不需要的頻率，從而大大影響音質，變調后會有部分失真[43]。另外，比較典型的方法是利用短時傅里葉變換原理，估計出短時幀的瞬時頻率，再乘以伸縮系數進行頻譜伸縮[44]。
參量法中最具代表性的方法是基於正弦模型原理。正弦模型[45]是由Quatier等人在1980年提出，它是目前應用最廣泛的語音模型。該模型將信號看作是一系列隨時間變化的正弦信號疊加。 很顯然，時間規整后瞬時頻率不變，保證了音調不變，但是時間過程擴展為原來的倍。很顯然，變調不變速處理后，各個頻率成分隨系數拉伸或者收縮。對應於濁音，為隨時間變化的第一諧波，即基頻；其他頻率成分對應於其它諧波。
由上分析可知，基於正弦模型的變調方法最大難點在於提高諧波分析的精確度，降低參數估計的復雜度[46]。
變速不變調方法
語音變速不變調，即語音時長規整，是指不改變原說話人的音調及語義信息，只改變說話人的語速。
語音變速不變調算法有三大類：時域法、頻域法、參量法，如表2-1所示。
表2-1 變速不變調算法分類
| 時域法 |
頻域法 |
參量法 |
| 剪貼法 |
LSEE-MSTFTM |
相位聲碼器 |
| SOLA、SOLA-FS |
|
正弦模型 |
| TD-PSOLA |
|
|
時域法包括：剪貼法、同步波形疊加法(Synchronized Overlap-Add, SOLA)、固定同步波形疊加法(Synchronized Overlap-Add and Fixed Synthesis, SOLAFS)、時域基音同步疊加法(Time-Domain Pitch Synchronized Overlap-Add, TD-PSOLA) 波形相似疊加法(waveform similarity overlap-and-add, WSOLA)等。
剪貼法[32]由Fairbanks等人在1958年提出，其基本思想是將語音划分為若干連續不重疊幀，然后重復或者丟棄其中某些語音幀，來實現語速變慢或者加快。
這種方法只是簡單的重復或者丟棄語音幀，這樣會造成相鄰兩幀之間波形不連續，基音發生斷裂，因此語音質量較差。
為了減小波形不連續現象，Gabor等人改善了剪貼法，在相鄰兩幀部分進行平滑處理，但是基音斷裂現象仍然存在。
為了在改善波形不連續的同時，減小基音斷裂現象，S.Roucos等人在提出了同步波形疊加法(SOLA)[33]。該算法主要分為分解、合成兩個階段。
分解階段完成原始語音信號的分幀任務，為了減小不連續現象，一般在分幀的同時進行加窗平滑處理，分解后的所有幀用於合成變速語音，這里假設幀長為N，分析幀移為Sa。合成階段則可以分成兩個步驟：
第一步，確定初步合成重疊距離。按照變速因子a=Ss/Sa，改變分解階段相鄰兩幀的幀移距離，Ss=Sa*a，即保持分解階段的第一幀位置不變，移動之后的各幀，使得相鄰兩幀的距離由為Sa變為Ss，這樣便可獲得初步合成幀。
第二步，確定最終合成幀的起始位置。如果將初步合成幀直接進行疊加合成，則會造成基音斷裂。為了減小該現象，通過在已合成第m幀第Ss個采樣點的鄰域[-kmax, kmax]內移動搜索和分解階段第m幀信號的波形相關最大的位置k(m)，如式2-1，確定最終合成幀的初始位置，保證疊加部分波形相似，減小基頻斷裂。
(2-1)
其中，k(m) [-kmax, kmax]，表示搜索偏移量，表示原始語音信號，表示合成變速后的輸出語音信號，L表示分解階段第m幀信號與已合成信號的重疊區域長度。
由於偏移量k(m)的存在，所以SOLA算法不能精確的改變語音的時間長度。
為了解決該問題，D.J.Hejna提出了固定同步波形疊加法(SOLAFS)[34]，該算法與SOLA算法很相似，不同的是在合成時固定了合成步長Ss，而在分解階段第m幀的鄰域[-kmax, kmax]內移動搜索與已合成第m幀信號波形相關最大的位置k(m)。即使式2-2最大。
為了進一步達到基音同步的目的，Moulines提出了基音同步波形疊加算法(PSOLA)[35]，著名的語音分析軟件Praat的變速功能就是基於此原理[36]。在該算法中，變速和變調是兩個獨立的過程，由不同的參數控制。其主要做法是首先進行基音同步分析，標記基音周期一系列位置點，以這些標記為中心，將原始語音划分成若干合成語音單元。通過重復或者省略合成單元來實現語速時長的控制，通過改變相鄰合成單元的重疊長度或者重采樣結合變速來改變語音的基頻。該算法主要缺點有兩點，首先基音周期及其起始點的判定誤差會嚴重影響PSOLA技術的效果，另外計算量大，很難滿足實時的變速與變調處理。
同樣地，為了減小基音斷裂和相位不連續問題，Verhelst和Roelands 提出了波形相似疊加法(WSOLA)[37]，該方法計算量低於PSOLA，同時輸出的語音質量高。Grofit對該方法進行了改進，使其也適用於音樂信號的變速處理[38]。
頻域法中最具代表性的方法是LSEE-MSTFTM(The Least-Square Error Estimation From the Modified Short-Time Fourier Transform Magnitude)[39]，該算法是基於短時傅里葉變換來實現的，利用最小均方誤差原則，尋找一個時域信號的短時傅里葉變換幅度譜逼近理想變速信號的頻譜。
參數法是指對語音信號建立模型，然后根據需要修改模型相關參數來達到改變音調和語速的目的。它包括相位聲碼器[40]和正弦模型法。相位聲碼器是將語音通過帶通濾波器組分解成若干正弦信號，然后對正弦信號隨時間變化的幅度和相位通過內插和抽取進行時域壓擴，最后經過合成便完成變速不變調。正弦模型法與相位聲碼器方法相似，需要估算出模型的瞬時幅度、瞬時頻率、瞬時相位等參量，合成效果較好，但是復雜度較高。
參考這個鏈接的一個例子。用的是變調不變速技術，主要是頻域相位伸縮+重采樣。
code:
function y = pvoc(x, r, n) % y = pvoc(x, r, n) Time-scale a signal to r times faster with phase vocoder % x is an input sound. n is the FFT size, defaults to 1024. % Calculate the 25%-overlapped STFT, squeeze it by a factor of r, % inverse spegram. % 2000-12-05, 2002-02-13 dpwe@ee.columbia.edu. Uses pvsample, stft, istft % $Header: /home/empire6/dpwe/public_html/resources/matlab/pvoc/RCS/pvoc.m,v 1.3 2011/02/08 21:08:39 dpwe Exp $ if nargin < 3 n = 1024; end % With hann windowing on both input and output, % we need 25% window overlap for smooth reconstruction hop = n/4; % Effect of hanns at both ends is a cumulated cos^2 window (for % r = 1 anyway); need to scale magnitudes by 2/3 for % identity input/output %scf = 2/3; % 2011-02-07: this factor is now included in istft.m scf = 1.0; % Calculate the basic STFT, magnitude scaled X = scf * stft(x', n, n, hop); % Calculate the new timebase samples [rows, cols] = size(X); t = 0:r:(cols-2); % Have to stay two cols off end because (a) counting from zero, and % (b) need col n AND col n+1 to interpolate % Generate the new spectrogram X2 = pvsample(X, t, hop); % Invert to a waveform y = istft(X2, n, n, hop)';
pvsample:phase vecoder sample
function c = pvsample(b, t, hop) % c = pvsample(b, t, hop) Interpolate an STFT array according to the 'phase vocoder' % b is an STFT array, of the form generated by 'specgram'. % t is a vector of (real) time-samples, which specifies a path through % the time-base defined by the columns of b. For each value of t, % the spectral magnitudes in the columns of b are interpolated, and % the phase difference between the successive columns of b is % calculated; a new column is created in the output array c that % preserves this per-step phase advance in each bin. % hop is the STFT hop size, defaults to N/2, where N is the FFT size % and b has N/2+1 rows. hop is needed to calculate the 'null' phase % advance expected in each bin. % Note: t is defined relative to a zero origin, so 0.1 is 90% of % the first column of b, plus 10% of the second. % 2000-12-05 dpwe@ee.columbia.edu % $Header: /homes/dpwe/public_html/resources/matlab/dtw/../RCS/pvsample.m,v 1.3 2003/04/09 03:17:10 dpwe Exp $ if nargin < 3 hop = 0; end [rows,cols] = size(b); N = 2*(rows-1); if hop == 0 % default value hop = N/2; end % Empty output array c = zeros(rows, length(t)); % Expected phase advance in each bin dphi = zeros(1,N/2+1); dphi(2:(1 + N/2)) = (2*pi*hop)./(N./(1:(N/2))); % Phase accumulator % Preset to phase of first frame for perfect reconstruction % in case of 1:1 time scaling ph = angle(b(:,1)); % Append a 'safety' column on to the end of b to avoid problems % taking *exactly* the last frame (i.e. 1*b(:,cols)+0*b(:,cols+1)) b = [b,zeros(rows,1)]; ocol = 1; for tt = t % Grab the two columns of b bcols = b(:,floor(tt)+[1 2]); tf = tt - floor(tt); bmag = (1-tf)*abs(bcols(:,1)) + tf*(abs(bcols(:,2))); % calculate phase advance dp = angle(bcols(:,2)) - angle(bcols(:,1)) - dphi'; % Reduce to -pi:pi range dp = dp - 2 * pi * round(dp/(2*pi)); % Save the column c(:,ocol) = bmag .* exp(j*ph); % Cumulate phase, ready for next frame ph = ph + dphi' + dp; ocol = ocol+1; end
STFT/ISTFT
function D = stft(x, f, w, h, sr)
% D = stft(X, F, W, H, SR) Short-time Fourier transform.
% Returns some frames of short-term Fourier transform of x. Each
% column of the result is one F-point fft (default 256); each
% successive frame is offset by H points (W/2) until X is exhausted.
% Data is hann-windowed at W pts (F), or rectangular if W=0, or
% with W if it is a vector.
% Without output arguments, will plot like sgram (SR will get
% axes right, defaults to 8000).
% See also 'istft.m'.
% dpwe 1994may05. Uses built-in 'fft'
% $Header: /home/empire6/dpwe/public_html/resources/matlab/pvoc/RCS/stft.m,v 1.4 2010/08/13 16:03:14 dpwe Exp $
if nargin < 2; f = 256; end
if nargin < 3; w = f; end
if nargin < 4; h = 0; end
if nargin < 5; sr = 8000; end
% expect x as a row
if size(x,1) > 1
x = x';
end
s = length(x);
if length(w) == 1
if w == 0
% special case: rectangular window
win = ones(1,f);
else
if rem(w, 2) == 0 % force window to be odd-len
w = w + 1;
end
halflen = (w-1)/2;
halff = f/2; % midpoint of win
halfwin = 0.5 * ( 1 + cos( pi * (0:halflen)/halflen));
win = zeros(1, f);
acthalflen = min(halff, halflen);
win((halff+1):(halff+acthalflen)) = halfwin(1:acthalflen);
win((halff+1):-1:(halff-acthalflen+2)) = halfwin(1:acthalflen);
end
else
win = w;
end
w = length(win);
% now can set default hop
if h == 0
h = floor(w/2);
end
c = 1;
% pre-allocate output array
d = zeros((1+f/2),1+fix((s-f)/h));
for b = 0:h:(s-f)
u = win.*x((b+1):(b+f));
t = fft(u);
d(:,c) = t(1:(1+f/2))';
c = c+1;
end;
% If no output arguments, plot a spectrogram
if nargout == 0
tt = [0:size(d,2)]*h/sr;
ff = [0:size(d,1)]*sr/f;
imagesc(tt,ff,20*log10(abs(d)));
axis('xy');
xlabel('time / sec');
ylabel('freq / Hz')
% leave output variable D undefined
else
% Otherwise, no plot, but return STFT
D = d;
end
%================
function x = istft(d, ftsize, w, h)
% X = istft(D, F, W, H) Inverse short-time Fourier transform.
% Performs overlap-add resynthesis from the short-time Fourier transform
% data in D. Each column of D is taken as the result of an F-point
% fft; each successive frame was offset by H points (default
% W/2, or F/2 if W==0). Data is hann-windowed at W pts, or
% W = 0 gives a rectangular window (default);
% W as a vector uses that as window.
% This version scales the output so the loop gain is 1.0 for
% either hann-win an-syn with 25% overlap, or hann-win on
% analysis and rect-win (W=0) on synthesis with 50% overlap.
% dpwe 1994may24. Uses built-in 'ifft' etc.
% $Header: /home/empire6/dpwe/public_html/resources/matlab/pvoc/RCS/istft.m,v 1.5 2010/08/12 20:39:42 dpwe Exp $
if nargin < 2; ftsize = 2*(size(d,1)-1); end
if nargin < 3; w = 0; end
if nargin < 4; h = 0; end % will become winlen/2 later
s = size(d);
if s(1) ~= (ftsize/2)+1
error('number of rows should be fftsize/2+1')
end
cols = s(2);
if length(w) == 1
if w == 0
% special case: rectangular window
win = ones(1,ftsize);
else
if rem(w, 2) == 0 % force window to be odd-len
w = w + 1;
end
halflen = (w-1)/2;
halff = ftsize/2;
halfwin = 0.5 * ( 1 + cos( pi * (0:halflen)/halflen));
win = zeros(1, ftsize);
acthalflen = min(halff, halflen);
win((halff+1):(halff+acthalflen)) = halfwin(1:acthalflen);
win((halff+1):-1:(halff-acthalflen+2)) = halfwin(1:acthalflen);
% 2009-01-06: Make stft-istft loop be identity for 25% hop
win = 2/3*win;
end
else
win = w;
end
w = length(win);
% now can set default hop
if h == 0
h = floor(w/2);
end
xlen = ftsize + (cols-1)*h;
x = zeros(1,xlen);
for b = 0:h:(h*(cols-1))
ft = d(:,1+b/h)';
ft = [ft, conj(ft([((ftsize/2)):-1:2]))];
px = real(ifft(ft));
x((b+1):(b+ftsize)) = x((b+1):(b+ftsize))+px.*win;
end;
測試code,因為變為Q = 5倍數，參數設置1/Q：
[d,sr]=wavread('aa.wav');
e = pvoc(d, 0.2);
f = resample(e,1,5); % NB: 0.2 = 1/5
結果
