>> % 2001-03-06-sinmodel.diary
>> % Examples of using sinewave modeling tools
>> % dpwe@ee.columbia.edu
>> 
>> % Load piano example
>> [d,r] = wavread('gould8.wav');
>> 
>> % Calculate spectrogram and 'instantaneous frequency'-gram (ifgram)
>> [I,S] = ifgram(d,256,256,128,r);
>> % FFT len 256 with 256 pt window, 128 point advance
>> 
>> % Track the local maxima in the spectrogram,
>> % (includes track-breaking at large magnitude increases)
>> [R,M] = extractrax(abs(S));
col 58 (l=54) tk 1: large step of 6.221191 at 27.341796, breaking track
col 58 (l=53) tk 5: large step of 4.426359 at 33.908660, breaking track
...
col 274 (l=55) tk 2: large step of 3.465624 at 15.630228, breaking track
col 326 (l=41) tk 4: large step of 6.605876 at 21.057121, breaking track
col 331 (l=84) tk 1: large step of 2.355255 at 40.604342, breaking track
>> % Each row of R and M describes a track.  Each column corresponds 
>> % to a time-frame column in spectrogram S.  For columns where track 
>> % is undefined, R value is NaN (which makes the point be skipped for 
>> % plotting).  Otherwise, R is the interpolated bin index of the 
>> % spectrogam peak, and M is the interpolated magnitude. 
>> size(R)

ans =

    60   335

>> % After pruning away very short tracks, there were 60 tracks 
>> % returned, extending over 335 time frames
>>
>> % To see the frequency tracks plotted on time-frequency axes:
>> plot(R')
>> % To have axes as time in seconds versus freq in Hz
>> tt = [0:334]*128/r;
>> % 128/r because each frame is a 128 sample advance, with r samples per sec
>> plot(tt, (R-1)*r/256)
>> % R-1 because bin indexes start at 1 to be Matlab-like
>> % r/256 because the 256 point FFT has the last bin correspond to r
>> 
>> % Can show superimposed on spectrogram:
>> specgram(d,256,r);
>> hold on
>> plot(tt, (R-1)*r/256)
>> % Clearly picked out most of the peaks
>> 
>> % However, for best reconstruction we should actually re-estimate 
>> % the track frequencies by interpolating the 'instantaneous frequency' 
>> % surface at each row position
>> F = colinterpvals(R,I);
>> % (colinterpvals does interpolation on each column of a surface)
>> % Similarly, if we want to *cancel out* the harmonic part to leave 
>> % just the untracked residual, we'll need to match the phase of 
>> % each track exactly, so interpolate the phases:
>> P = -colinterpvals(R, unwrap(angle(S)));
>> % have to negate the phase because resynthesis wants phase relative 
>> % to t=0, but analysis returns phase relative to analysis frame
>> % unwrap mostly not needed, just in case adjacent bins straddle wrap point
>> % (unwrap acts along columns, as we want).
>> 
>> % Pad frq, mag and phs matrices with an extra column at each end to 
>> % account for the half-windows lost by having NFFT=2*HOP
>> F = [F(:,1)*0,F,F(:,335)*0];
>> % Multiplying adjacent column by zero preserves NaNs for inactive rows
>> M = [M(:,1)*0,M,M(:,335)*0];
>> P = [P(:,1)*0,P,P(:,335)*0];
>> % Resynthesize sinusoids for each track.  Have to scale M by 2 because 
>> % we only looked at the positive frequency half
>> sws = synthphtrax(F,M*2,P,r,128);
>> size(sws)

ans =

           1       43008

>> size(d)

ans =

       43058           1

>> % Calculate residual
>> dr = d(1:43008) - sws(1:43008)';
>> subplot(311)
>> specgram(d,256,r);
>> subplot(312)
>> specgram(sws,256,r);
>> subplot(313)
>> specgram(dr,256,r);
>> soundsc(d,r)
>> soundsc(sws,r)
>> soundsc(dr,r)
>> % Residual is mainly noisy onset transients as notes hit, plus 
>> % a few low-energy harmonics that didn't get tracked
>> 
>> % Could now search through tracks to find groups e.g. by common onset.
>> % Individual tracks are picked out by selecting rows from F, M
>> subplot(111)
>> plot(M(1:10,10))
>> % At time frame 10, the track in row 2 has the highest magnitude
>> % Let's look at it in detail:
>> plot(M(2,:))
>> % Only exists for first 60 bins or so
>> % Also, magnitude easier to read in dB
>> ii = 1:65;
>> plot3(tt(ii),F(2,ii),20*log10(M(2,ii)))
>> grid
>> % Frequency jumps about a bit suspiciously, but only over about 0.5%
>> % Magnitude decay pretty linear
>> % Onset values definitely weird
>>
>> % Plot all the early tracks in 3-D
>> plot3(tt(ii),F([1:8],ii),20*log10(M([1:8],ii)))
>> grid
>> view([-30 60])
>> % Gives you a picture of several simultaneous harmonics
>> % Add projections onto t-f plane
>> hold on
>> plot3(tt(ii),F([1:8],ii),-80*ones(1,65))
>> 
>> diary