diff extra/sinemodelpeakint.m @ 10:6840f77b83aa

commit

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/extra/sinemodelpeakint.m	Sun Apr 21 10:55:35 2013 +0200
@@ -0,0 +1,78 @@
+function [y]=sinemodelpeakint(x,w,N,t)
+%initializing values
+M = length(w); % window size - the longer the more frequency resolution
+N2 = N/2+1; % positive part of the spectrum
+Ns= 1024; % FFT size for synthesis (even)
+H = 256; % analysis/synthesishop size
+soundlength = length(x); % length of input sound array - samples
+
+fftbuffer = zeros(N,1); % initialize buffer for FFT
+
+%Create a loop to step through the sound array x
+%initializing the loop
+hNs = Ns/2; % half synthesis window size
+hM = (M-1)/2; % half analysis window size used to overlap windows
+
+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
+
+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);
+
+while pin<pend
+    xw = x(pin-hM:pin+hM).*w(1:M)'; % window the input sound - STFT definition
+  
+    %zero phased window
+    fftbuffer(:) = 0; % reset buffer
+    fftbuffer(1:(M+1)/2) = xw((M+1)/2:M); % zero-phase fftbuffer
+    fftbuffer(N-(M-1)/2+1:N) = xw(1:(M-1)/2);
+
+    %compute FFT of the zero phased frame
+    X = fft(fftbuffer);
+
+    %calculate  magnitude and phase spectrum of of positive frequencies
+    mX = 20*log10(abs(X(1:N2)));
+    pX = unwrap(angle(X(1:N2)));
+    
+        
+    %Find the locations, ploc, of the local maxima above a given 
+    %threshold, t, in each magnitude spectrum by finding changes of slope.
+    ploc = 1+find((mX(2:N2-1)>t).*(mX(2:N2-1)>mX(3:N2)).*(mX(2:N2-1)>mX(1:N2-2)));
+    
+    %Find the magnitudes, pmag, and phases, pphase, of the obtained
+    %locations.
+    pmag = mX(ploc);
+    %pmag = mX(ploc)*0.4; 
+    pphase = pX(ploc); 
+    
+    
+    %peak interpolation
+    [iploc, ipmag, ipphase] = peakinterp (mX, pX, ploc);
+    
+    
+% plotting the peak values on top of the magnitude and phase spectra at the
+% magnitude and phase frame
+  
+subplot(211); 
+plot(mX);
+hold on
+plot(ploc,pmag,'m*');   %m is for magenta color
+plot(nploc,npmag,'c*'); %plot the new magnitude spectra
+title('magnitude and peak values');
+hold off
+subplot(212);
+plot(pX);
+hold on
+plot(ploc,pphase,'m*');
+plot(nploc,npphase,'c*');
+title('phase and peak values');
+hold off
+pause
+ 
\ No newline at end of file