|
yading@10
|
1 function y = sinemodel(x, w, N, t)
|
|
yading@10
|
2 % Analysis/synthesis of a sound using the sinusoidal model
|
|
yading@10
|
3 % x: input sound, w: analysis window (odd size), N: FFT size,
|
|
yading@10
|
4 % t: threshold in negative dB, y: output sound
|
|
yading@10
|
5 M = length(w); % analysis window size
|
|
yading@10
|
6 Ns= 1024; % FFT size for synthesis (even)
|
|
yading@10
|
7 H = 256; % analysis/synthesishop size
|
|
yading@10
|
8 N2= N/2+1; % size of positive spectrum
|
|
yading@10
|
9 soundlength = length(x); % length of input sound array
|
|
yading@10
|
10 hNs = Ns/2; % half synthesis window size
|
|
yading@10
|
11 hM = (M-1)/2; % half analysis window size
|
|
yading@10
|
12 pin = max(N2,1+hM); % initialize sound pointer to middle of analysis window
|
|
yading@10
|
13 pend = soundlength-max(N2,hM); % last sample to start a frame
|
|
yading@10
|
14 fftbuffer = zeros(N,1); % initialize buffer for FFT
|
|
yading@10
|
15 y = zeros(soundlength,1); % initialize output array
|
|
yading@10
|
16 w = w/sum(w); % normalize analysis window
|
|
yading@10
|
17 sw = zeros(Ns,1);
|
|
yading@10
|
18 ow = triang(2*H-1); % overlapping window (triangular window to avoid too much overlapping)
|
|
yading@10
|
19 ovidx = Ns/2+1-H+1:Ns/2+H; % overlap indexes
|
|
yading@10
|
20 sw(ovidx) = ow(1:2*H-1);
|
|
yading@10
|
21 bh = blackmanharris(Ns); % synthesis window
|
|
yading@10
|
22 bh = bh ./ sum(bh); % normalize synthesis window
|
|
yading@10
|
23 sw(ovidx) = sw(ovidx) ./ bh(ovidx);
|
|
yading@10
|
24 while pin<pend
|
|
yading@10
|
25 %-----analysis-----%
|
|
yading@10
|
26 xw = x(pin-hM:pin+hM).*w(1:M); % window the input sound
|
|
yading@10
|
27 %zero phased window
|
|
yading@10
|
28 fftbuffer(:)=0;
|
|
yading@10
|
29 fftbuffer(1:(M+1)/2) = xw((M+1)/2:M); % zero-phase window in fftbuffer
|
|
yading@10
|
30 fftbuffer(N-(M-1)/2+1:N) = xw(1:(M-1)/2);
|
|
yading@10
|
31
|
|
yading@10
|
32 X = fft(fftbuffer); % compute the FFT
|
|
yading@10
|
33 mX = 20*log10(abs(X(1:N2))); % magnitude spectrum of positive frequencies
|
|
yading@10
|
34 pX = unwrap(angle(X(1:N/2+1))); % unwrapped phase spectrum
|
|
yading@10
|
35 ploc = 1 + find((mX(2:N2-1)>t) .* (mX(2:N2-1)>mX(3:N2)).* (mX(2:N2-1)>mX(1:N2-2))); % find peaks
|
|
yading@10
|
36 [ploc,pmag,pphase] = peakinterp(mX,pX,ploc); % refine peak values
|
|
yading@10
|
37 %-----synthesis-----%
|
|
yading@10
|
38 plocs = (ploc-1)*Ns/N+1; % adapt peak locations to synthesis FFT
|
|
yading@10
|
39 Y = genspecsines(plocs,pmag,pphase,Ns); % generate spec sines
|
|
yading@10
|
40 yw = fftshift(real(ifft(Y))); % time domain of sinusoids
|
|
yading@10
|
41 y(pin-hNs:pin+hNs-1) = y(pin-hNs:pin+hNs-1) + sw.*yw(1:Ns); % overlap-add
|
|
yading@10
|
42 pin = pin+H; % advance the sound pointer
|
|
yading@10
|
43 end |
