boblsturm@0: function [x,CF,CM] = chromsynth(F,bp,sr)
boblsturm@0: % [x,CF,CM] = chromsynth(F,bp,sr)
boblsturm@0: %   Resynthesize a chroma feature vector to audio
boblsturm@0: %   F is 12 rows x some number of columns, one per beat
boblsturm@0: %   bp is the period of one beat in sec (or a vector of beat times)
boblsturm@0: %   sr is the sampling rate of the output waveform
boblsturm@0: %   x is returned as a 12 semitone-spaced sines modulated by F
boblsturm@0: %   CF,CM return actual sinusoid matrices passed to synthtrax
boblsturm@0: %   Actual Shepard tones now implemented! 2007-04-19
boblsturm@0: % 2006-07-14 dpwe@ee.columbia.edu
boblsturm@0: 
boblsturm@0: if nargin < 2; bp = 0.5; end % 120 bpm
boblsturm@0: if nargin < 3; sr = 22050; end
boblsturm@0: 
boblsturm@0: [nchr,nbts] = size(F);
boblsturm@0: 
boblsturm@0: % resynth
boblsturm@0: if length(bp) == 1
boblsturm@0:   bups = 8; % upsampling factor
boblsturm@0:   framerate = bups/bp;
boblsturm@0:   ncols = nbts*bups;
boblsturm@0:   CMbu = zeros(nchr, ncols);
boblsturm@0:   for i = 1:bups
boblsturm@0:     CMbu = CMbu + upsample(F', bups, i-1)';
boblsturm@0:   end
boblsturm@0: else
boblsturm@0:   % vector of beat times - quantize
boblsturm@0:   framerate = 50; % frames per sec
boblsturm@0:   nbeats = length(bp);
boblsturm@0:   lastbeat = bp(nbeats) + (bp(nbeats) - bp(nbeats-1));
boblsturm@0:   ncols = round(lastbeat * framerate);
boblsturm@0:   CMbu = zeros(nchr, ncols);
boblsturm@0:   xF = [zeros(12,1),F];
boblsturm@0:   for i = 1:ncols
boblsturm@0:     CMbu(:,i) = xF(:,max(find(i/framerate >= [0,bp])));
boblsturm@0:   end
boblsturm@0: end
boblsturm@0:   
boblsturm@0: %CFbu = repmat(440*2.^([0:(nchr-1)]'/nchr),1,ncols);
boblsturm@0: %x = synthtrax(CFbu,CMbu,sr,round(sr/framerate));
boblsturm@0: octs = 7;
boblsturm@0: basefrq = 27.5;  % A1; +6 octaves = 3520
boblsturm@0: CFbu = repmat(basefrq*2.^([0:(nchr-1)]'/nchr),1,ncols);
boblsturm@0: CF = [];
boblsturm@0: CM = [];
boblsturm@0: % what bin is the center freq?
boblsturm@0: f_ctr = 440;
boblsturm@0: f_sd = 0.5;
boblsturm@0: f_bins = basefrq*2.^([0:(nchr*octs - 1)]/nchr);
boblsturm@0: f_dist = log(f_bins/f_ctr)/log(2)/f_sd;  
boblsturm@0: % actually just = ([0:(nchr*octs - 1)]/nchr-log2(f_ctr/basefrq))/f_sd
boblsturm@0: % Gaussian weighting centered of f_ctr, with f_sd
boblsturm@0: f_wts = exp(-0.5*f_dist.^2);
boblsturm@0: for oct = 1:octs
boblsturm@0:   CF = [CF;(2^oct)*CFbu];
boblsturm@0:   % pick out particular weights
boblsturm@0:   CM = [CM;diag(f_wts((oct-1)*nchr+[1:nchr]))*CMbu];
boblsturm@0: end
boblsturm@0: % Remove sines above nyquist
boblsturm@0: CFok = (CF(:,1) < sr/2);
boblsturm@0: CF = CF(CFok,:);
boblsturm@0: CM = CM(CFok,:);
boblsturm@0: % Synth the sines
boblsturm@0: x = synthtrax(CF,CM,sr,round(sr/framerate));
boblsturm@0: 
boblsturm@0: % Playing synth along with audio:
boblsturm@0: %>> rdc = chromsynth(Dy.F(:,160+[1:300]),Dy.bts(160+[1:300]) - Dy.bts(160),sr);
boblsturm@0: %>> Dy.bts(161)*sr              
boblsturm@0: %ans =
boblsturm@0: %      279104
boblsturm@0: %>> size(rdc)
boblsturm@0: %ans =
boblsturm@0: %           1      498401
boblsturm@0: %>> ddd = db(279104+[1:498401]);
boblsturm@0: %>> soundsc([ddd,rdc'/40],sr)