音频变调技术
今天看到群里有人讨论这个问题，记录一下。
主要内容转自：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
结果
