Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- function y = harmonicmodel(x, w, N, t, fs)
- % Authors: J. Bonada, X. Serra, X. Amatriain, A. Loscos
- % Analysis/synthesis of a sound using the sinusoidal model
- % x: input sound, w: analysis window (odd size), N: FFT size,
- % t: threshold in negative dB, y: output sound
- %
- %--------------------------------------------------------------------------
- % This source code is provided without any warranties as published in
- % DAFX book 2nd edition, copyright Wiley & Sons 2011, available at
- % http://www.dafx.de. It may be used for educational purposes and not
- % for commercial applications without further permission.
- %--------------------------------------------------------------------------
- M = length(w); % analysis window size
- Ns= 1024; % FFT size for synthesis (even)
- H = 256; % analysis/synthesishop size
- N2= N/2+1; % size of positive spectrum
- soundlength = length(x); % length of input sound array
- hNs = Ns/2; % half synthesis window size
- hM = (M-1)/2; % half analysis window size
- pin = max(H+1,1+hM); % initialize sound pointer to middle of analysis window
- pend = soundlength-max(H,hM); % last sample to start a frame
- fftbuffer = zeros(N,1); % initialize buffer for FFT
- y = zeros(soundlength,1); % initialize output array
- w = w/sum(w); % normalize analysis window
- sw = zeros(Ns,1);
- ow = triang(2*H-1); % overlapping window
- ovidx = Ns/2+1-H+1:Ns/2+H; % overlap indexes
- sw(ovidx) = ow(1:2*H-1);
- bh = blackmanharris(Ns); % synthesis window
- bh = bh ./ sum(bh); % normalize synthesis window
- sw(ovidx) = sw(ovidx) ./ bh(ovidx);
- k = [0:1:N];
- freq = fs*k/N;
- minf0 = 100;
- maxf0 = 20000;
- ef0max = 100;
- while pin<pend
- %-----analysis-----%
- xw = x(pin-hM:pin+hM).*w(1:M); % window the input sound
- fftbuffer(:) = 0; % reset buffer
- fftbuffer(1:(M+1)/2) = xw((M+1)/2:M); % zero-phase window in fftbuffer
- fftbuffer(N-(M-1)/2+1:N) = xw(1:(M-1)/2);
- X = fft(fftbuffer); % compute the FFT
- mX = 20*log10(abs(X(1:N2))); % magnitude spectrum of positive frequencies
- pX = unwrap(angle(X(1:N/2+1))); % unwrapped phase spectrum
- 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
- [ploc,pmag,pphase] = peakinterp(mX,pX,ploc); % refine peak values
- % f0 = f0detection(mX, fs, ploc, pmag, ef0max, minf0, maxf0);
- f0 = f0detectionyin(x, fs, 1200, minf0, maxf0) ;
- nH=20;
- maxhd=.2;
- hf = (f0>0).*(f0.*(1:nH)); % initialize harmonic frequencies
- hi = 1; % initialize harmonic index
- npeaks = length(ploc); % number of peaks found
- while (f0>0 && hi<=nH && hf(hi)<fs/2) % find harmonic peaks
- [dev,pei] = min(abs((ploc(1:npeaks)-1)/N*fs-hf(hi))); % closest peak
- if ((hi==1 || ~any(hloc(1:hi-1)==ploc(pei))) && dev<maxhd*hf(hi))
- hloc(hi) = ploc(pei); % harmonic locations
- hmag(hi) = pmag(pei); % harmonic magnitudes
- hphase(hi) = pphase(pei); % harmonic phases
- end
- hi = hi+1; %increase harmonic index
- end
- % f0
- % f0yin
- % subplot(2,1,1);
- % plot(freq(1:N2),mX);
- % hold on;
- % plot(f0,mX(round(f0*N/fs)),'g*');
- % plot(fs*hloc/N,hmag,'rX');
- % plot(f0yin,mX(round(f0yin*N/fs)),'g*');
- % hold off;
- % pause;
- % subplot(2,1,2);
- % plot(f0,mX(round(f0*N/fs)),'X');
- %-----synthesis-----%
- plocs = (ploc-1)*Ns/N; % adapt peak locations to synthesis FFT
- Y = genspecsines(plocs,pmag,pphase,Ns); % generate spec sines
- yw = fftshift(real(ifft(Y))); % time domain of sinusoids
- %(pin-hNs:pin+hNs-1)
- y(pin-hNs:pin+hNs-1) = y(pin-hNs:pin+hNs-1) + sw.*yw(1:Ns); % overlap-add
- pin = pin+H; % advance the sound pointer
- end
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement