Mercurial > hg > pmhd
changeset 4:aeccfa83c7a4
Merge
| author | Katerina <katkost@gmail.com> |
|---|---|
| date | Sat, 20 Apr 2013 13:03:57 +0100 |
| parents | 1c0f36c348d4 (diff) e0a7176da80a (current diff) |
| children | 840a4c453edc |
| files | |
| diffstat | 40 files changed, 1472 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra codes/D.m Sat Apr 20 13:03:57 2013 +0100 @@ -0,0 +1,5 @@ +function y = D(x,N) +% Calculate rectangular window transform (Dirichlet kernel) +y = sin(N*x/2)./sin(x/2); +y(find(y~=y))=N; % avoid NaN if x==0 +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra codes/FFT of a Simple Sinusoid.m Sat Apr 20 13:03:57 2013 +0100 @@ -0,0 +1,45 @@ +% Example 1: FFT of a DFT_Sinusoid + +% Parameters: +N=64; % Must be a power of two +T=1; % Set sampling rate to 1 +A=1; % Sinusoidal Amplitude +phi=0; % Sinusoidal phase +f=0.25; % Frequency (cycles/sample) +n=[0:N-1]; % Discrete Time axis +x=A*cos(2*pi*n*f*T+phi); % Sampled sinusoid +X=fft(x); % Spectrum + +% Plot time data: +figure(1); +subplot(3,1,1); +plot(n,x,'*k'); +ni= [0:.1:N-1]; % Interpolated time axis +hold on; +plot(ni,A*cos(2*pi*ni*f*T+phi),'-k'); grid off; +title('Sinusoid at 1/4 the Spampling Rate'); +xlabel('Time (samples)'); +ylabel('Amplitude'); +text(-8,1,'a)'); +hold off; + +% Plot spectral magnitude; +magX= abs(X); +fn=[0:1/N:1-1/N]; %Normalized frequency axis +sublplot(3,1,2); +stem(fn,magX,'ok');grid on; +xlabel('Normalized Frequency(cycles per sample))'); +ylabel('Magnitude(Linear)'); +text(-.11,40,'b)'); + +% Same thing on a dB scale: +spec=20*log10(magX); %Spectral magitude in dB +subplot(3,1,3); +plot(fn,spec,'--ok'); grid on; +axis([0 1 -350 50]); +xlabel('Normalized Frequency (ccles per sample))'); +ylabel('Magnitude(dB)'); +text('Magnitude(dB)'); +text(-.11,50, 'c)'); +cmd=['print -deps ', '../eps/example1.eps']; +disp(cmd); eval(cmd); \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra codes/Lab2 1b.m Sat Apr 20 13:03:57 2013 +0100 @@ -0,0 +1,1 @@ +x2=[0 1 0 0 0 0 0 0]; \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra codes/Lab2 1d.m Sat Apr 20 13:03:57 2013 +0100 @@ -0,0 +1,4 @@ +x4=[2 1 0 0 0 0 0 1 ]; +X4=fft(x4,8); +ejex=linspace(0,1,8); +plot(ejex,abs(X4),'*'); \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra codes/Lab21b.m Sat Apr 20 13:03:57 2013 +0100 @@ -0,0 +1,12 @@ +x3=[1 1 1 1 1 1 1 1 ]; +X3=fft(x3,8); +ejex=linspace(0,1,8); +figure('Name','X3'); +subplot(131); +plot(ejex,abs(X3),'*'); +subplot(132); +plot(ejex,angle(X3),'*'); + + + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra codes/Lab21d.m Sat Apr 20 13:03:57 2013 +0100 @@ -0,0 +1,4 @@ +x2=[0 1 0 0 0 0 0 0]; +X2=fft(x2,8); +plot(x2); +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra codes/Lab4 stptpeaks.m Sat Apr 20 13:03:57 2013 +0100 @@ -0,0 +1,51 @@ +function y=stptpeaks(x, w, N, H, t) +% Analysis/synthesis of a sound using the peaks +% of the short-time fourier transform +% x: input sound, +% w: analysis window (odd size), +% N: FFT size, +% H: hop size, +% t: threshold in negative dB, +% y: output sound + +M = length(w); % analysis window size +N2 = N/2+1; % size of positive spectrum +soundlength = length(x); % length of input sound array +hM = (M-1)/2; % half analysis window size + +pin = 1+hM; % initialize sound pointer at the middle of analysis window +pend = soundlength-hM; % last sample to start a frame + +fftbuffer = zeros(N,1); % initialize buffer for FFT +yw = zeros(M,1); % initialize output sound frame +y = zeros(soundlength,1); % initialize output array + +w = w/sum(w); % normalize analysis window +sw = hanning(M); % synthesis window +sw = sw./sum(sw); + +while pin<pend +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 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:N2))); % 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))); % peakss +pmag = mX(ploc); % magnitude of peaks +pphase = pX(ploc); % phase of peaks +num_peak=length(ploc)/2; %Just positive peaks + +axis=linspace(0,pi,N2); %Axis for the whole frequency plot +axis_peak=axis(ploc(1:num_peak));%Axis for just the peaks plot + +subplot(2,1,1); plot(axis,mX); hold on; %Magnitude plot +plot(axis_peak,pmag(1:num_peak),'rx');%Magnitude peaks +subplot(2,1,2); plot(axis,pX); hold on; %Phase plot +plot(axis_peak,pphase(1:num_peak),'rx');%Phase peaks + +pin = pin+H; % advance sound pointer +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra codes/Lab4.m Sat Apr 20 13:03:57 2013 +0100 @@ -0,0 +1,44 @@ +function y = stpt(x, w, N, H, t) +% Analysis/synthesis of a sound using the peaks +% of the short-time fourier transform +% x: input sound, +% w: analysis window (odd size), +% N: FFT size, +% H: hop size, +% t: threshold in negative dB, +% y: output sound +M = length(w); % analysis window size +N2 = N/2+1; % size of positive spectrum +soundlength = length(x); % length of input sound array +hM = (M-1)/2; % half analysis window size +pin = 1+hM; % initialize sound pointer at the middle of analysis window +pend = soundlength-hM; % last sample to start a frame +fftbuffer = zeros(N,1); % initialize buffer for FFT +yw = zeros(M,1); % initialize output sound frame +y = zeros(soundlength,1); % initialize output array +w = w/sum(w); % normalize analysis window +sw = hanning(M); % synthesis window +sw = sw./sum(sw); +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 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:N2))); % 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))); % peakss +pmag = mX(ploc); % magnitude of peaks +pphase = pX(ploc); % phase of peaks +%-----synthesis-----% +Y = zeros(N,1); % initialize output spectrum +Y(ploc) = 10.^(pmag/20).*exp(i.*pphase); % generate positive freq. +Y(N+2-ploc) = 10.^(pmag/20).*exp(-i.*pphase); % generate negative freq. +fftbuffer = real(ifft(Y)); % real part of the inverse FFT +yw((M+1)/2:M) = fftbuffer(1:(M+1)/2); % undo zero phase window +yw(1:(M-1)/2) = fftbuffer(N-(M-1)/2+1:N); +y(pin-hM:pin+hM) = y(pin-hM:pin+hM) + H*N*sw.*yw(1:M); % overlap-add +pin = pin+H; % advance sound pointer +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra codes/TWM.m Sat Apr 20 13:03:57 2013 +0100 @@ -0,0 +1,55 @@ + +function [f0, f0error] = TWM (ploc, pmag, N, fs, minf0, maxf0) +% Two-way mismatch algorithm (by Beauchamp&Maher) +% ploc, pmag: peak locations and magnitudes, N: size of complex spectrum +% fs: sampling rate of sound, minf0: minimum f0, maxf0: maximum f0 +% f0: fundamental frequency detected, f0error: error measure +pfreq = (ploc-1)/N*fs; % frequency in Hertz of peaks +[zvalue,zindex] = min(pfreq); +if (zvalue==0) % avoid zero frequency peak + pfreq(zindex) = 1; + pmag(zindex) = -100; +end +ival2 = pmag; +[Mmag1,Mloc1] = max(ival2); % find peak with maximum magnitude +ival2(Mloc1) = -100; % clear max peak +[Mmag2,Mloc2]= max(ival2); % find second maximum magnitude peak +ival2(Mloc2) = -100; % clear second max peak +[Mmag3,Mloc3]= max(ival2); % find third maximum magnitude peak +nCand = 3; % number of possible f0 candidates for each max peak +f0c = zeros(1,3*nCand); % initialize array of candidates +f0c(1:nCand)=(pfreq(Mloc1)*ones(1,nCand))./((nCand+1-(1:nCand))); % candidates +f0c(nCand+1:nCand*2)=(pfreq(Mloc2)*ones(1,nCand))./((nCand+1-(1:nCand))); +f0c(nCand*2+1:nCand*3)=(pfreq(Mloc3)*ones(1,nCand))./((nCand+1-(1:nCand))); +f0c = f0c((f0c<maxf0)&(f0c>minf0)); % candidates within boundaries +if (isempty(f0c)) % if no candidates exit + f0 = 0; f0error=100; + return +end +harmonic = f0c; +ErrorPM = zeros(fliplr(size(harmonic))); % initialize PM errors +MaxNPM = min(10,length(ploc)); +for i=1:MaxNPM % predicted to measured mismatch error + difmatrixPM = harmonic' * ones(size(pfreq))'; + difmatrixPM = abs(difmatrixPM-ones(fliplr(size(harmonic)))*pfreq'); + [FreqDistance,peakloc] = min(difmatrixPM,[],2); + Ponddif = FreqDistance .* (harmonic'.^(-0.5)); + PeakMag = pmag(peakloc); + MagFactor = 10.^((PeakMag-Mmag1)./20); + ErrorPM = ErrorPM+(Ponddif+MagFactor.*(1.4*Ponddif-0.5)); + harmonic = harmonic+f0c; +end +ErrorMP = zeros(fliplr(size(harmonic))); % initialize MP errors +MaxNMP = min(10,length(pfreq)); +for i=1:length(f0c) % measured to predicted mismatch error + nharm = round(pfreq(1:MaxNMP)/f0c(i)); + nharm = (nharm>=1).*nharm + (nharm<1); + FreqDistance = abs(pfreq(1:MaxNMP) - nharm*f0c(i)); + Ponddif = FreqDistance.* (pfreq(1:MaxNMP).^(-0.5)); + PeakMag = pmag(1:MaxNMP); + MagFactor = 10.^((PeakMag-Mmag1)./20); + ErrorMP(i) = sum(MagFactor.*(Ponddif+MagFactor.*(1.4*Ponddif-0.5))); +end +Error = (ErrorPM/MaxNPM) + (0.3*ErrorMP/MaxNMP); % total errors +[f0error, f0index] = min(Error); % get the smallest error +f0 = f0c(f0index); % f0 with the smallest error \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra codes/ejempolo de tres matrics HeBeR.m Sat Apr 20 13:03:57 2013 +0100 @@ -0,0 +1,5 @@ +x3=[1 1 1 1 1 1 1 1 ]; +X3=fft(x3,8); +ejex=linspace(0,1,8); +plot(ejex,abs(X3),'*'); +%plot(ejex,angle(X3),'*'); \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra codes/f0detection.m Sat Apr 20 13:03:57 2013 +0100 @@ -0,0 +1,15 @@ +function f0 = f0detection(mX, fs, ploc, pmag, ef0max, minf0, maxf0) +% Fundamental frequency detection function +% mX: magnitude spectrum, fs: sampling rate, ploc, pmag: peak loc and mag, +% ef0max: maximim error allowed, minf0: minimum f0, maxf0: maximum f0 +% f0: fundamental frequency detected +N = length(mX)*2; % size of complex spectrum +nPeaks = length(ploc); % number of peaks +f0 = 0; % initialize output +if(nPeaks>3) % at least 3 peaks in spectrum for trying to find f0 + nf0peaks = min(50,nPeaks); % use a maximum of 50 peaks + [f0,f0error] = TWM(ploc(1:nf0peaks),pmag(1:nf0peaks),N,fs,minf0,maxf0); + if (f0>0 && f0error>ef0max) % limit the possible error by ethreshold + f0 = 0; + end +end;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra codes/f0detectionyin.m Sat Apr 20 13:03:57 2013 +0100 @@ -0,0 +1,60 @@ +function f0 = f0detectionyin(x,fs,ws,minf0,maxf0) +% x: input signal, fs: sampling rate, ws: integration window length, minf0: minimum f0, maxf0: maximum f0 +% f0: fundamental frequency detected in Hz +maxlag = ws-2; % maximum lag +th = 0.1; % set threshold +d = zeros(maxlag,1); % init variable +d2 = zeros(maxlag,1); % init variable +% compute d(tau) +x1 = x(1:ws); +cumsumx = sum(x1.^2); +cumsumxl = cumsumx; +xy = xcorr(x(1:ws*2),x1); +xy = xy(ws*2+1:ws*3-2); +for lag=0:maxlag-1 + d(1+lag) = cumsumx + cumsumxl - 2*xy(1+lag); + cumsumxl = cumsumxl - x(1+lag).^2 + x(1+lag+ws+1)^2; +end +cumsum = 0; +% compute d'(tau) +d2(1) = 1; +for lag=1:maxlag-1 + cumsum = cumsum + d(1+lag); + d2(1+lag) = d(1+lag)*lag./cumsum; +end +% limit the search to the target range +minf0lag = 1+round(fs./minf0); % compute lag corresponding to minf0 +maxf0lag = 1+round(fs./maxf0); % compute lag corresponding to maxf0 +if (maxf0lag>1 && maxf0lag<maxlag) + d2(1:maxf0lag) = 100; % avoid lags shorter than maxf0lag +end +if (minf0lag>1 && minf0lag<maxlag) + d2(minf0lag:end) = 100; % avoid lags larger than minf0lag +end +% find the best candidate +mloc = 1 + find((d2(2:end-1)<d2(3:end)).*(d2(2:end-1)<d2(1:end-2))); % minima +candf0lag = 0; +if (length(mloc)>0) + I = find(d2(mloc)<th); + if (length(I)>0) + candf0lag = mloc(I(1)); + else + [Y,I2] = min(d2(mloc)); + candf0lag = mloc(I2); + end + candf0lag = candf0lag; % this is zero-based indexing + if (candf0lag>1 & candf0lag<maxlag) + % parabolic interpolation + lval = d2(candf0lag-1); + val = d2(candf0lag); + rval= d2(candf0lag+1); + candf0lag = candf0lag + .5*(lval-rval)./(lval-2*val+rval); + end +end +ac = min(d2); +f0lag = candf0lag-1; % convert to zero-based indexing +f0 = fs./f0lag; % compute candidate frequency in Hz +if (ac > 0.2) % voiced/unvoiced threshold + f0 = 0; % set to unvoiced +end + \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra codes/genbh92lobe.m Sat Apr 20 13:03:57 2013 +0100 @@ -0,0 +1,18 @@ +function y = genbh92lobe(x) +% Calculate transform of the Blackman-Harris 92dB window +% x: bin positions to compute (real values), y: transform values +N = 512; +f = x*pi*2/N; % frequency sampling +df = 2*pi/N; +y = zeros(size(x)); % initialize window +consts = [.35875, .48829, .14128, .01168]; % window constants +for m=0:3 + y = y + consts(m+1)/2*(D(f-df*m,N)+D(f+df*m,N)); % sum Dirichlet kernels +end +y = y/N/consts(1); % normalize +end +function y = D(x,N) +% Calculate rectangular window transform (Dirichlet kernel) +y = sin(N*x/2)./sin(x/2); +y(find(y~=y))=N; % avoid NaN if x==0 +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra codes/genspecsines.m Sat Apr 20 13:03:57 2013 +0100 @@ -0,0 +1,26 @@ +function Y = genspecsines (ploc, pmag, pphase, N) +% Compute a spectrum from a series of sine values +% iploc, ipmag, ipphase: sine locations, magnitudes and phases, N: size of complex +% spectrum; Y: generated complex spectrum of sines +Y =zeros(N,1); % initialize output spectrum +hN = N/2+1; % size of positive freq. spectrum +for i=1:length(ploc); % generate all sine spectral lobes + loc = ploc(i)-1; % location of peak (zero-based indexing),range ]0,hN-1[ + if (loc<=1||loc>=hN-1) continue; end; % avoid frequencies out of range + binremainder = round(loc)-loc; + lb = [binremainder-4:binremainder+4]'; % main lobe (real value) bins to read + lmag = genbh92lobe(lb)*10.^(pmag(i)/20); % lobe magnitudes of the complex exponential + b = 1+[round(loc)-4:round(loc)+4]'; % spectrum bins to fill (1-based indexing) + for m=1:9 + if (b(m)<1) % peak lobe crosses DC bin + Y(2-b(m)) = Y(2-b(m)) + lmag(m)*exp(-1i*pphase(i)); + elseif (b(m)>hN) % peak lobe croses Nyquist bin + Y(2*hN-b(m)) = Y(2*hN-b(m)) + lmag(m)*exp(-1i*pphase(i)); + else % peak lobe in positive freq. range + Y(b(m)) = Y(b(m)) + lmag(m)*exp(1i*pphase(i)) ... + + lmag(m)*exp(-1i*pphase(i))*(b(m)==1||b(m)==hN); + end + end + Y(hN+1:end) = conj(Y(hN-1:-1:2)); % fill the rest of the spectrum +end +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra codes/genspecsines2.m Sat Apr 20 13:03:57 2013 +0100 @@ -0,0 +1,11 @@ +Y = genspecsines([10.3, 20.6, 30.9], [-10, -13, -16], [0, 0, 0], 1024); +%Y_sin = genspecsines_sin([10.3, 20.6, 30.9], [-10, -13, -16],[0, 0, 0], 1024); +figure; +subplot(2,2,1) +plot(Y); +subplot(2,2,2) +plot(Y_sin); +subplot(2,2,3) +plot(ifft(Y)); +subplot(2,2,4) +plot(ifft(Y_sin)); \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra codes/genspecsines3.m Sat Apr 20 13:03:57 2013 +0100 @@ -0,0 +1,12 @@ +Y = genspecsines([10.3, 20.6, 30.9], [-10, -13, -16], [0, 0, 0], 1024); +%plot(ifft(Y)); +f1=10.3; +f2=20.6; +f3=30.9; +t = [0:1/44100:0.1-1/44100]; +x1 = cos(2*pi*f1*t); +x2 = cos(2*pi*f2*t); +x3 = cos(2*pi*f3*t); +x = x1 + x2+ x3; +plot(x); +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra codes/genspecsines_sin.m Sat Apr 20 13:03:57 2013 +0100 @@ -0,0 +1,14 @@ +function Y = genspecsines_sin(ploc, pmag, pphase, N) +y=zeros(N,1); +w=blackmanharris(N); +w=w./sum(w); +w(2:N)=w(2:N)*2; +for k=1:length(ploc) +f=ploc(k)/N; +A=10ˆ(pmag(k)/20); +phi=pphase(k); +n=[0:1:N-1]'; +y=A*cos(2*pi*f.*n+phi)+y; +end +y=y.*w; +Y=abs(fft(y)); \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra codes/gong.m Sat Apr 20 13:03:57 2013 +0100 @@ -0,0 +1,41 @@ +function y = stpt(x, w, N, H, t) +% Analysis/synthesis of a sound using the peaks +% of the short-time fourier transform +% x: input sound, w: analysis window (odd size), N: FFT size, H: hop size, +% t: threshold in negative dB, y: output sound +M = length(w); % analysis window size +N2 = N/2+1; % size of positive spectrum +soundlength = length(x); % length of input sound array +hM = (M-1)/2; % half analysis window size +pin = 1+hM; % initialize sound pointer at the middle of analysis window +pend = soundlength-hM; % last sample to start a frame +fftbuffer = zeros(N,1); % initialize buffer for FFT +yw = zeros(M,1); % initialize output sound frame +y = zeros(soundlength,1); % initialize output array +w = w/sum(w); % normalize analysis window +sw = hanning(M); % synthesis window +sw = sw./sum(sw); +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 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:N2))); % 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))); % peaks +pmag = mX(ploc); % magnitude of peaks +pphase = pX(ploc); % phase of peaks +%-----synthesis-----% +Y = zeros(N,1); % initialize output spectrum +Y(ploc) = 10.^(pmag/20).*exp(i.*pphase); % generate positive freq. +Y(N+2-ploc) = 10.^(pmag/20).*exp(-i.*pphase); % generate negative freq. +fftbuffer = real(ifft(Y)); % real part of the inverse FFT +yw((M+1)/2:M) = fftbuffer(1:(M+1)/2); % undo zero phase window +yw(1:(M-1)/2) = fftbuffer(N-(M-1)/2+1:N); +y(pin-hM:pin+hM) = y(pin-hM:pin+hM) + H*N*sw.*yw(1:M); % overlap-add +pin = pin+H; % advance sound pointer +end +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra codes/gong_notes.m Sat Apr 20 13:03:57 2013 +0100 @@ -0,0 +1,14 @@ +Spectral Peaks + +GONG: + +sound(y,16000) +y=spectralpeaks(gong, hanning(101) , 1020, 51, -80) %pierde armonicos +y=spectralpeaks(gong, hanning(501) , 1020, 51, -80) %suena casi igual que el original +y=spectralpeaks(gong, hanning(501) , 1020, 251, -80) %it is even smoother + +y=spectralpeaks(gong, hanning(2001), 2001, 501, -80) %sounds good +y=spectralpeaks(gong, hanning(2001), 2001, 501, -40) %excite some harmonic components + + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra codes/harmonicmodel.m Sat Apr 20 13:03:57 2013 +0100 @@ -0,0 +1,62 @@ +function y = harmonicmodel(x, fs, w, N, t, nH, minf0, maxf0, f0et, maxhd) +% Analysis/synthesis of a sound using the sinusoidal harmonic model +% x: input sound, fs: sampling rate, w: analysis window (odd size), +% N: FFT size (minimum 512), t: threshold in negative dB, +% nH: maximum number of harmonics, minf0: minimum f0 frequency in Hz, +% maxf0: maximim f0 frequency in Hz, +% f0et: error threshold in the f0 detection (ex: 5), % maxhd: max. relative deviation in harmonic detection (ex: .2) +% y: output sound +M = length(w); % analysis window size +Ns= 1024; % FFT size for synthesis +H = 256; % hop size for analysis and synthesis +N2 = N/2+1; % size postive 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(hNs+1,1+hM); % initialize sound pointer to middle of analysis window +pend = soundlength-hM; % last sample to start a frame +fftbuffer = zeros(N,1); % initialize buffer for FFT +y = zeros(soundlength+Ns/2,1); % output sound +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 + %-----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 + 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,f0et,minf0,maxf0); % find f0 + hloc = zeros(nH,1); % initialize harmonic locations + hmag = zeros(nH,1)-100; % initialize harmonic magnitudes + hphase = zeros(nH,1); % initialize harmonic phases + 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 + hloc(1:hi-1) = (hloc(1:hi-1)~=0).*((hloc(1:hi-1)-1)*Ns/N+1); % synth. locs + %-----synthesis-----% + Yh = genspecsines(hloc(1:hi-1),hmag,hphase,Ns); % generate sines + yh = fftshift(real(ifft(Yh))); % sines in time domain + y(pin-hNs:pin+hNs-1) = y(pin-hNs:pin+hNs-1) + sw.*yh(1:Ns); % overlap-add + pin = pin+H; % advance the input sound pointer +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra codes/harmonicmodel1.m Sat Apr 20 13:03:57 2013 +0100 @@ -0,0 +1,90 @@ +function [y]=harmonicmodel1(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(hNs+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-hNs+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); + + + %plot for a window + if (pin==1+hM+20*H) + figure + subplot(2,1,1) + plot(mX) + hold on + plot(iploc,mX(ploc),'*') + title('magnitude and peak values'); + hold off + subplot(2,1,2) + plot(pX) + hold on + plot(iploc,pX(ploc),'*') + title('phase and peak values'); + hold off + + %number of peaks + npeaksm = length(pmag) + npeaksp = length(pphase) + end + + + %-----synthesis-----% +plocs = (ploc-1)*Ns/N+1; % adapt peak locations to synthesis FFT +Y = genspecsines(plocs,pmag,pphase,Ns); % generate spec sines +yw = fftshift(real(ifft(Y))); % time domain of sinusoids +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 +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra codes/lab2_4.m Sat Apr 20 13:03:57 2013 +0100 @@ -0,0 +1,9 @@ +N=1024; +t = [0:1/44100:0.1-1/44100]; +x1=cos(2*pi*100*44100/1024*t); +x801 = x1(1:801); +xlength=length(x801); +fftbuffer=zeros(N, 1); +fftbuffer(1:(xlength+1)/2) = x801((xlength+1)/2:xlength); +fftbuffer(N-(xlength-1)/2+1:N) = x801(1: (xlength-1)/2); +plot(unwrap(angle(fftbuffer,1024))); \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra codes/lab5.m Sat Apr 20 13:03:57 2013 +0100 @@ -0,0 +1,63 @@ +clear all +%1.1 +%% Piano sound + +% hfile = '4192__RealRhodesSounds__D4.wav'; +% [x3, Fs3, nbits3, readinfo3] = wavread(hfile); +% %sound(x3, Fs); +% y = sinemodel_(x3, hamming(513)', 1024, 100,-35); +% wavwrite(y, Fs3,'__RealRhodesSounds_D4.wav'); + + +%% 2.1 Compute the magnitude values on the to of the main lobe +% loc=0; +% mag=0; +% +% binremainder = round(loc)-loc; +% %Blackman-Harris 92dB +% lb = [binremainder-4:binremainder+4]'; % main lobe (real value) bins to read +% lmag = genbh92lobe(lb,1024)*10.^(mag/20); % lobe magnitudes +% figure +% plot(lmag); + + +%% 2.2 Compute a spectrum of size N from a series of sinusoidal values + +% Y = genspecsines([10.3, 20.6, 30.9], [-10, -13, -16], [0, 0, 0], 1024); +% +% %inverse fft +% y=fftshift(real(ifft(Y))); +% plot(y); +% +% %the sinusoids +% Y1 = (-10)*exp(1*i);y1=ifft(Y1); +% Y2 = (-13)*exp(1*i);y2=fftshift(real(ifft(Y2))); +% Y3 = (-16)*exp(1*i);y3=fftshift(real(ifft(Y3))); +% %add the sinusoids in the time domain +% ytime = y1+y2+y3; +% figure +% plot(y1); + +%% Fugue + +% hfile = '22125__Connum__memento_to_Bachs_fugue_II.wav'; +% [x3, Fs3, nbits3, readinfo3] = wavread(hfile); +% %sound(x3, Fs); +% y = sinemodel(x3, hamming(2001)', 1024, -60); +% wavwrite(y, Fs3,'fugue1.wav'); + +% %% Soprano +% +% hfile = '30084__HerbertBoland__FemalePhrase1.wav'; +% [x3, Fs3, nbits3, readinfo3] = wavread(hfile); +% %sound(x3, Fs); +% y = sinemodel(x3, hamming(2001), 1024, -50); +% wavwrite(y, Fs3,'female1.wav'); + +%% Soprano + +hfile = '32848__Quilt__Quilty_s_Old_School_Piano.wav'; +[x3, Fs3, nbits3, readinfo3] = wavread(hfile); +%sound(x3, Fs); +y = sinemodel_(x3, hamming(2001)', 1024, -50); +wavwrite(y, Fs3,'piano1.wav');
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra codes/peakinterp.m Sat Apr 20 13:03:57 2013 +0100 @@ -0,0 +1,11 @@ +function [iploc, ipmag, ipphase] = peakinterp(mX, pX, ploc) +% Parabolic interpolation of spectral peaks +% mX: magnitude spectrum, pX: phase spectrum, ploc: locations of peaks +% iploc, ipmag, ipphase: interpolated values +% note that ploc values are assumed to be between 2 and length(mX)-1 +val = mX(ploc); % magnitude of peak bin +lval = mX(ploc-1); % magnitude of bin at left +rval= mX(ploc+1); % magnitude of bin at right +iploc = ploc + .5*(lval-rval)./(lval-2*val+rval); % center of parabola +ipmag = val-.25*(lval-rval).*(iploc-ploc); % magnitude of peaks +ipphase = interp1(1:length(pX),pX,iploc,'linear'); % phase of peaks \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra codes/peaks.m Sat Apr 20 13:03:57 2013 +0100 @@ -0,0 +1,51 @@ +%Hecho con Álvaro +function y=peaks(x,w,N,H,t) +%Analysis/sysntesis of a sound using the peaks +%of the short-time fourier transform +%x: input sound +%w: analysis window (odd size) +%N: FFT size +%H: hop size +%t: treshold in negative dB +%y: output sound + +%Internal variables +M=length(w); +N2=N/2+1; +soundlength=length(x); +hM=(M-1)/2; + +%Pointers +pin=1+hM; +pend=soundlength-hM; + +%Initialize buffers +fftbuffer=zeros(N,1); +yw=zeros(M,1); +y=zeros(soundlength,1); + +while pin<pend + xw=x(pin-hM:pin+hM).*w(1:M); %Frame selector + fftbuffer=zeros(N,1); %Zero-phase windowing + fftbuffer(1:(M+1)/2)=xw((M+1)/2:M); + fftbuffer(N-(M-1)/2+1:N)=xw(1:(M-1)/2); + + X=fft(fftbuffer); %FFT Performance + mX=20*log10(abs(X(1:N2))); %Compute FFT Magnitude + pX=phase(X(1:N2)); %Compute FFT Phase + + ploc=1+find((mX(2:N2-1)>t).*(mX(2:N2-1)>mX(3:N2)).*(mX(2:N2-1)>mX(1:N2-2))); %Peak search (levels above threshold) + pmag=mX(ploc); %Find the peak magnitude + pphase=pX(ploc); %Find the peak phase + % posipeak=length(ploc)/2; %Just positive peaks + + axis=linspace(0,pi,abs(N2)); %Axis for the whole frequency plotting + axis_peak=axis(ploc);%Axis for just the peaks plotting + + subplot(2,1,1); plot(axis,mX); hold on ; %Magnitude plotting + plot(axis_peak,pmag,'rx');hold off;%Magnitude peaks plot red x’s + subplot(2,1,2); plot(axis,pX); hold on ;%Phase plotting + plot(axis_peak,pphase,'rx');hold off;%Phase peaks plot red x’s + getframe; + pin=pin+H; %Advance the pointer +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra codes/peakseaker.m Sat Apr 20 13:03:57 2013 +0100 @@ -0,0 +1,50 @@ +function y=peaks(x,w,N,H,t) +%Analysis/sysntesis of a sound using the peaks +%of the short-time fourier transform +%x: input sound +%w: analysis window (odd size) +%N: FFT size +%H: hop size +%t: treshold in negative dB +%y: output sound + +%Internal variables +M=length(w); +N2=N/2+1; +soundlength=length(x); +hM=(M-1)/2; + +%Pointers +pin=1+hM; +pend=soundlength-hM; + +%Initialize buffers +fftbuffer=zeros(N,1); +yw=zeros(M,1); +y=zeros(soundlength,1); + +while pin<pend + xw=x(pin-hM:pin+hM).*w(1:M); %Frame selector + fftbuffer=zeros(N,1); %Zero-phase windowing + fftbuffer(1:(M+1)/2)=xw((M+1)/2:M); + fftbuffer(N-(M-1)/2+1:N)=xw(1:(M-1)/2); + + X=fft(fftbuffer); %FFT Performance + mX=20*log10(abs(X(1:N2))); %Compute FFT Magnitude + pX=phase(X(1:N2)); %Compute FFT Phase + + ploc=1+find((mX(2:N2-1)>t).*(mX(2:N2-1)>mX(3:N2)).*(mX(2:N2-1)>mX(1:N2-2))); %Peak search (levels above threshold) + pmag=mX(ploc); %Find the peak magnitude + pphase=pX(ploc); %Find the peak phase + posipeak=length(ploc)/2; %Just positive peaks + + axis=linspace(0,pi,abs(N2)); %Axis for the whole frequency plotting + axis_peak=axis(ploc(1:num_peak));%Axis for just the peaks plotting + + subplot(2,1,1); plot(axis,mX); hold off ; %Magnitude plotting + plot(axis_peak,pmag(1:posipeak),'rx');%Magnitude peaks plot red x’s + subplot(2,1,2); plot(axis,pX); hold off ;%Phase plotting + plot(axis_peak,pphase(1:posipeak),'rx');%Phase peaks plot red x’s + + pin=pin+H; %Advance the pointer +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra codes/plot_peaks.m Sat Apr 20 13:03:57 2013 +0100 @@ -0,0 +1,14 @@ +function plot_peaks(Xmag,Xphase,ploc,t) +figure(1); +subplot(2,1,1) +plot(Xmag); +hold on; +plot(ploc,Xmag(ploc),'r*'); +plot([1,length(Xmag)],[t,t],'g'); +hold off; +subplot(2,1,2); +plot(Xphase); +hold on; +plot(ploc,Xphase(ploc),'r*'); +hold off; +pause; \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra codes/sinemodel.m Sat Apr 20 13:03:57 2013 +0100 @@ -0,0 +1,43 @@ +function y = sinemodel(x, w, N, t) +% 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 +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(N2,1+hM); % initialize sound pointer to middle of analysis window +pend = soundlength-max(N2,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 (triangular window to avoid too much overlapping) +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 +%-----analysis-----% +xw = x(pin-hM:pin+hM).*w(1:M); % window the input sound +%zero phased window +fftbuffer(:)=0; +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 +%-----synthesis-----% +plocs = (ploc-1)*Ns/N+1; % adapt peak locations to synthesis FFT +Y = genspecsines(plocs,pmag,pphase,Ns); % generate spec sines +yw = fftshift(real(ifft(Y))); % time domain of sinusoids +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 \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra codes/sinemodel1.m Sat Apr 20 13:03:57 2013 +0100 @@ -0,0 +1,90 @@ +function [y]=sinemodel1(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); + + + %plot for a window + if (pin==1+hM+20*H) + figure + subplot(2,1,1) + plot(mX) + hold on + plot(iploc,mX(ploc),'*') + title('magnitude and peak values'); + hold off + subplot(2,1,2) + plot(pX) + hold on + plot(iploc,pX(ploc),'*') + title('phase and peak values'); + hold off + + %number of peaks + npeaksm = length(pmag) + npeaksp = length(pphase) + end + + + %-----synthesis-----% +plocs = (ploc-1)*Ns/N+1; % adapt peak locations to synthesis FFT +Y = genspecsines(plocs,pmag,pphase,Ns); % generate spec sines +yw = fftshift(real(ifft(Y))); % time domain of sinusoids +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 +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra codes/sinemodel_.m Sat Apr 20 13:03:57 2013 +0100 @@ -0,0 +1,79 @@ + +function [y]=sinemodel_(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= 2048; % FFT size for synthesis (even) +H = 512; % 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))); % unwrapped phase spectrum + + + %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))); %peaks + + %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); + + %plot for a window + + subplot(2,1,1) + plot(mX) + hold on + plot(ploc,pmag,'*'); + plot(iploc,ipmag,'c'); + title('magnitude peak values'); + hold off + subplot(2,1,2) + plot(pX) + hold on + plot(ploc,pphase,'m'); + plot(iploc,ipphase,'c*') + title('phase peak values'); + hold off + + %number of peaks + +end +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra codes/sinemodelpeakint.m Sat Apr 20 13:03:57 2013 +0100 @@ -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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra codes/spectralpeaks.m Sat Apr 20 13:03:57 2013 +0100 @@ -0,0 +1,41 @@ +function y = stpt(x, w, N, H, t) +% Analysis/synthesis of a sound using the peaks +% of the short-time fourier transform +% x: input sound, w: analysis window (odd size), N: FFT size, H: hop size, +% t: threshold in negative dB, y: output sound +M = length(w); % analysis window size +N2 = N/2+1; % size of positive spectrum +soundlength = length(x); % length of input sound array +hM = (M-1)/2; % half analysis window size +pin = 1+hM; % initialize sound pointer at the middle of analysis window +pend = soundlength-hM; % last sample to start a frame +fftbuffer = zeros(N,1); % initialize buffer for FFT +yw = zeros(M,1); % initialize output sound frame +y = zeros(soundlength,1); % initialize output array +w = w/sum(w); % normalize analysis window +sw = hanning(M); % synthesis window +sw = sw./sum(sw); +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 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:N2))); % 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))); % peaks +pmag = mX(ploc); % magnitude of peaks +pphase = pX(ploc); % phase of peaks +%-----synthesis-----% +Y = zeros(N,1); % initialize output spectrum +Y(ploc) = 10.^(pmag/20).*exp(i.*pphase); % generate positive freq. +Y(N+2-ploc) = 10.^(pmag/20).*exp(-i.*pphase); % generate negative freq. +fftbuffer = real(ifft(Y)); % real part of the inverse FFT +yw((M+1)/2:M) = fftbuffer(1:(M+1)/2); % undo zero phase window +yw(1:(M-1)/2) = fftbuffer(N-(M-1)/2+1:N); +y(pin-hM:pin+hM) = y(pin-hM:pin+hM) + H*N*sw.*yw(1:M); % overlap-add +pin = pin+H; % advance sound pointer +end +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra codes/stft.m Sat Apr 20 13:03:57 2013 +0100 @@ -0,0 +1,43 @@ +function [mY, pY] = stft(x, w, N, H) +% Analysis/synthesis of a sound using the short-time fourier transform +% x: input sound, w: analysis window (odd size), N: FFT size, H: hop size +% y: output sound +M = length(w); % analysis window size +N2 = N/2+1; % size of positive spectrum +soundlength = length(x); % length of input sound array +hM = (M-1)/2; % half analysis window size +pin = 1+hM; % initialize sound pointer in middle of analysis window +pend = soundlength-hM; % last sample to start a frame +fftbuffer = zeros(N,1); % initialize buffer for FFT or reset buffer +w = w/sum(w); % normalize analysis window +nframes=mod(soundlength-hM,H); +frame=1; +%fftbuffer=zeros(N,1); +mY=zeros(N2,nframes); +pY=zeros(N2,nframes); + +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 FFT +mX = 20*log10(abs(X(1:N2))); % magnitude spectrum of positive frequencies +pX = unwrap(angle(X(1:N2))); % unwrapped phase spect. of positive freq. +mY(:,frame)=mX; +pY(:,frame)=pX; +frame=frame+1; +pin = pin+H; % advance sound pointer +end +% Copied from Carles +mY_norm=mY-(mY./mY).*max(max(mY)); +%Spectrogram mode respresentation +surface(1:frame-1, linspace(0,pi,N2),mY_norm(:,1:frame-1),'LineStyle','none'); +%Axes adjustment +xlim([1 frame-1]); +ylim([0 pi]); + + + +[x,fs,nbit]=wavread('Ocarina.wav');
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra codes/stft2.m Sat Apr 20 13:03:57 2013 +0100 @@ -0,0 +1,34 @@ +function y = stft2(x, w, N, H) +% Analysis/synthesis of a sound using the short-time fourier transform +% x: input sound, w: analysis window (odd size), N: FFT size, H: hop size +% y: output sound +M = length(w); % analysis window size +N2 = N/2+1; % size of positive spectrum +soundlength = length(x); % length of input sound array +hM = (M-1)/2; % half analysis window size +pin = 1+hM; % initialize sound pointer in middle of analysis window +pend = soundlength-hM; % last sample to start a frame +fftbuffer = zeros(N,1); % initialize buffer for FFT +yw = zeros(M,1); % initialize output sound frame +y = zeros(soundlength,1); % initialize output array +w = w/sum(w); % normalize analysis window +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 FFT +mX = 20*log10(abs(X(1:N2))); % magnitude spectrum of positive frequencies +pX = unwrap(angle(X(1:N2))); % unwrapped phase spect. of positive freq. +%-----synthesis-----% +Y = zeros(N,1); % initialize output spectrum +Y(1:N2) = 10.^(mX/20).*exp(i.*pX); % generate positive freq. +Y(N2+1:N) = 10.^(mX(N2-1:-1:2)/20).*exp(-i.*pX(N2-1:-1:2)); +% generate neg.freq. +fftbuffer = real(ifft(Y)); % inverse FFT +yw(1:(M-1)/2) = fftbuffer(N-(M-1)/2+1:N); % undo zero-phase window +yw((M+1)/2:M) = fftbuffer(1:(M+1)/2); +y(pin-hM:pin+hM) = y(pin-hM:pin+hM) + H*yw(1:M); % overlap-add +pin = pin+H; % advance sound pointer +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra codes/stft_peak.m Sat Apr 20 13:03:57 2013 +0100 @@ -0,0 +1,51 @@ +function X = stft_peak(x, w, N, H, threshold) +% Analysis/synthesis of a sound using the short-time fourier transform +% x: input sound, +% w: analysis window (odd number of samples required) +% N: FFT size, +% H: hop size +% y: output sound +lengthx = length(x); % length of sound signal +w = w/sum(w); % normalize window +W = length(w); % window size (odd) +W2 = (W-1)/2; % half window +N2 = N/2+1; % half FFT size + +p_start = 1+W2; % first frame pointer +p_end = lengthx-W2; % last frame pointer +fft_buffer = zeros(N,1); % FFT buffer +for g=p_start:H:p_end % advance hop size each loop + % ANALYSIS + xw = w(1:W)'.*x(g-W2:g+W2); % multiply window by input frame + fft_buffer(1:W2+1) = xw(W2+1:W); % adjust buffer for zero-phase computation + fft_buffer(N-W2+1:N) = xw(1:W2); + X = fft(fft_buffer); % FFT of current buffer + mX = 20*log10(abs(X(1:N2))); % magnitude spectrum + pX = unwrap(angle(X(1:N2))); % phase spectrum (unwrapped) + ploc = 1+find((mX(2:N2-1)>threshold).*(mX(2:N2-1)>mX(3:N2)).*(mX(2:N2-1)>mX(1:N2-2))); + pmag = mX(ploc); + pphase = pX(ploc); + plength = length(ploc); + if(mod((g-p_start)/H+1,(floor(p_end/H/5)))==0) % select only 5 equally distributed frames of the sound + figure + subplot(2,1,1) + plot(mX) % plot magnitude and phase spectrums + hold on + plot(ploc(1:plength),pmag(1:plength),'xr'); % add peaks to the plot + xlim([1 N2]) + ylim([-80 0]) + title('Magnitude spectrum of frame') + ylabel('Magnitude spectrum (dB)') + xlabel('Frequency (rad)') + display(['Peak positions: ' num2str(ploc')]); + hold off + subplot(2,1,2) + plot(pX) % same with the phase + hold on + plot(ploc(1:plength),pphase(1:plength),'xr'); + xlim([1 N2]) +% ylim([-pi pi]) +% pause() + hold off + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra codes/stft_primer_ejemplo.m Sat Apr 20 13:03:57 2013 +0100 @@ -0,0 +1,34 @@ +function y = stft(x, w, N, H) +% Analysis/synthesis of a sound using the short-time fourier transform +% x: input sound, w: analysis window (odd size), N: FFT size, H: hop size +% y: output sound +M = length(w); % analysis window size +N2 = N/2+1; % size of positive spectrum +soundlength = length(x); % length of input sound array +hM = (M-1)/2; % half analysis window size +pin = 1+hM; % initialize sound pointer in middle of analysis window +pend = soundlength-hM; % last sample to start a frame +fftbuffer = zeros(N,1); % initialize buffer for FFT +yw = zeros(M,1); % initialize output sound frame +y = zeros(soundlength,1); % initialize output array +w = w/sum(w); % normalize analysis window +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 FFT +mX = 20*log10(abs(X(1:N2))); % magnitude spectrum of positive frequencies +pX = unwrap(angle(X(1:N2))); % unwrapped phase spect. of positive freq. +%-----synthesis-----% +Y = zeros(N,1); % initialize output spectrum +Y(1:N2) = 10.^(mX/20).*exp(i.*pX); % generate positive freq. +Y(N2+1:N) = 10.^(mX(N2-1:-1:2)/20).*exp(-i.*pX(N2-1:-1:2)); +% generate neg.freq. +fftbuffer = real(ifft(Y)); % inverse FFT +yw(1:(M-1)/2) = fftbuffer(N-(M-1)/2+1:N); % undo zero-phase window +yw((M+1)/2:M) = fftbuffer(1:(M+1)/2); +y(pin-hM:pin+hM) = y(pin-hM:pin+hM) + H*yw(1:M); % overlap-add +pin = pin+H; % advance sound pointer +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra codes/stftlab4.m Sat Apr 20 13:03:57 2013 +0100 @@ -0,0 +1,51 @@ +function y = stftlab4(x, w, N, H, t) +% Analysis/synthesis of a sound using the peaks +% of the short-time fourier transform +% x: input sound, w: analysis window (odd size), N: FFT size, H: hop size, +% t: threshold in negative dB, y: output sound +M = length(w); % analysis window size +N2 = N/2+1; % size of positive spectrum +soundlength = length(x); % length of input sound array +hM = (M-1)/2; % half analysis window size +pin = 1+hM; % initialize sound pointer at the middle of analysis window +pend = soundlength-hM; % last sample to start a frame +fftbuffer = zeros(N,1); % initialize buffer for FFT +yw = zeros(M,1); % initialize output sound frame +y = zeros(soundlength,1); % initialize output array +w = w/sum(w); % normalize analysis window +sw = hanning(M); % synthesis window +sw = sw./sum(sw); +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 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:N2))); % 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))); % peaks +[nploc,npmag,npphase] = peakinterp(mX,pX,ploc); % refine peak values. i use new parameters +%so that i can plot the new values at the figures existed +pmag = mX(ploc); % finding the magnitude of the obtained locations (peaks) +pphase = pX(ploc); % finding the phase of the obtained locations (peaks) + +% 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 + +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra codes/stpt.m Sat Apr 20 13:03:57 2013 +0100 @@ -0,0 +1,39 @@ +function y= stpt(x,w,N,H,t) +M=length(w); % Window size +pin=0; % Start pointer +pend=length(x)-M; % End pointer +y=zeros(1,length(x)); % Inicial. output +npeaks=0; % Number of peaks +w=w./sum(w); % Analy. Windw. Norm. +w_out=hanning(M)'; % Synthesis window. +w_out=w_out./sum(w_out); % Synth. Windw. Norm. +while(pin<pend) % Until the end... + %---Analysis + x_w=x(pin+1:pin+M).*w(1:M); % Windowing signal + fftbuffer=[x_w, zeros(1,N-M)]; % Zero-padding + fftbuffer=circshift(fftbuffer,[0,-(M-1)/2]);% Zero-phase shift + X=fft(fftbuffer,N); % FFT + Xmag=20*log10(abs(X(1:N/2+1))); % Magnitude in dBs + Xphase=unwrap(angle(X(1:N/2+1))); % Phase + %---Transformation Xmag and Xphase + ploc = 1+find((Xmag(2:N/2)>t).*(Xmag(2:N/2)>Xmag(3:N/2+1)).*(Xmag(2:N/2)>Xmag(1:N/2-1))); % Find the local max. + npeaks=npeaks+length(ploc); + pphase=Xphase(ploc); % Save phase values + pmag=Xmag(ploc); % Save mag. values + %plot_peaks(Xmag,Xphase,ploc,t); % Plotting peaks + Xmag=zeros(1,N/2+1)-1000; % Modifying Xmag... + Xphase=zeros(1,N/2+1); % Modifying Xphase... + %ploc=round(ploc.*0.96); % Transposition (fx) + %ploc=10:10:10*length(ploc); % Harmonics mapping(fx) + Xmag(ploc)=pmag; % Add the mag maxs. + Xphase(ploc)=pphase; % Add the phase maxs. + %---Synthesis + Y(1:N/2+1) = 10.^(Xmag/20).*exp(1i.*Xphase);% Generate posit. freq. + Y(N/2+2:N) = conj(fliplr(Y(2:N/2))); % Generate negat. freq + fftbuffer_out=real(ifft(Y)); % Clean the signal + fftbuffer_out=circshift(fftbuffer_out,[0,+(M-1)/2]); + % Undo zero-phasing + y(pin+1:pin+M)=y(pin+1:pin+M)+H*N*fftbuffer_out(1:M).*w_out; + % Overlap-add method + pin=pin+H; % Advance pointer +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra codes/stpt2.m Sat Apr 20 13:03:57 2013 +0100 @@ -0,0 +1,53 @@ +function y=stpt2(x,w,N,H,t) +%Analysis/synthesis of a sound using the peaks +%of the short-time Fourier transform +%x: input sound, w: analysis window (odd size), N: FFT size, H: hop +%size, t: threshold in negative dB, y: output sound + +%Internal variables +M=length(w); +N2=N/2+1; +soundlength=length(x); +hM=(M-1)/2; + +%Pointers +pin=1+hM; +pend=soundlength-hM; + +%Initialize buffers +fftbuffer=zeros(N,1); +yw=zeros(M,1); +y=zeros(soundlength,1); + +%Analysis/Synthesis windows +sw=hanning(M); +sw=sw./sum(sw); +w=w./sum(w); + +while pin<pend + %Analysis + xw=x(pin-hM:pin+hM).*w(1:M); %Frame selector + fftbuffer=zeros(N,1); %Zero-phase windowing + fftbuffer(1:(M+1)/2)=xw((M+1)/2:M); + fftbuffer(N-(M-1)/2+1:N)=xw(1:(M-1)/2); + + X=fft(fftbuffer); %FFT performance + mX=20*log10(abs(X(1:N2))); %Compute the FFT Magnitude (dB) + pX=angle(X(1:N2)); %Compute the FFT phase + + ploc=1+find((mX(2:N2-1)>t).*(mX(2:N2-1)>mX(3:N2)).*(mX(2:N2-1)>mX(1:N2-2))); %Peak finder above threshold (t) + pmag=mX(ploc); %Peak magnitude recovery + pphase=pX(ploc); %Peak phase recovery + + %Synthesis + Y=zeros(N,1); %Memory reserve + Y(ploc)=10.^(pmag/2).*exp(1i.*pphase); %Spectrum reconstruction from peak info + Y(N+2-ploc)=10.^(pmag/20).*exp(-1i.*pphase); %Spectrum reconstruction negative freq. + fftbuffer=real(ifft(Y)); %Compute inverse FFT + yw((M+1)/2:M)=fftbuffer(1:(M+1)/2); %Zero-phase windowing undo + yw(1:(M-1)/2)=fftbuffer(N-(M-1)/2+1:N); + y(pin-hM:pin+hM)=y(pin-hM:pin+hM)+H*N*sw.*yw(1:M); %overlap-add + pin=pin+H; %Pointer advance +end + + \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra codes/stptpeaks.m Sat Apr 20 13:03:57 2013 +0100 @@ -0,0 +1,49 @@ +function y=stptpeaks(x, w, N, H, t) +% Analysis/synthesis of a sound using the peaks +% of the short-time fourier transform +% x: input sound, +% w: analysis window (odd size), +% N: FFT size, +% H: hop size, +% t: threshold in negative dB, +% y: output sound + +M = length(w); % analysis window size +N2=N/2+1; % size of positive spectrum +soundlength = length(x); % length of input sound array +hM = (M-1)/2; % half analysis window size + +pin = 1+hM; % initialize sound pointer at the middle of analysis window +pend = soundlength-hM; % last sample to start a frame + +fftbuffer = zeros(N,1); % initialize buffer for FFT +yw = zeros(M,1); % initialize output sound frame +y = zeros(soundlength,1); % initialize output array + + + +while pin<pend +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 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:N2))); % 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))); % peakss +pmag = mX(ploc); % magnitude of peaks +pphase = pX(ploc); % phase of peaks +posipeaks=length(ploc)/2; %Just positive peaks + +axis=linspace(0,pi,N2); %Axis for the whole frequency plot +axispeaks=axis(ploc(1:posipeaks));%Axis for just the peaks plot + +subplot(2,1,1); plot(axis,mX); hold off; %Magnitude plot +plot(axispeaks,pmag(1:posipeaks),'rx');%Magnitude peaks +subplot(2,1,2); plot(axis,pX); hold off; %Phase plot +plot(axispeaks,pphase(1:posipeaks),'rx');%Phase peaks + +pin = pin+H; % advance sound pointer +end \ No newline at end of file