I proivided those files for you in my answer to your previous question Unrecognized function or variable 'chromagram_IF'.
It likely would have been better if you had continued this in that thread rather than starting a new thresd.
If you provide your ‘audio_2024-05-12_17-59-56.ogg’ file, I might be able to jelp you with this. Use the ‘paperclip’ icon to upload it. If ‘.ogg’ is an unrecognised file type, then first ues the zip function to ‘zip’ the file, and then upload the .zip file.
EDIT — (13 May 2024 at 12:09)
You may need to use the chromsynth function. See the provided documentation under ‘Chroma Synthesis’.
EDIT — (13 May 2024 at 15:21)
Try this —
% load('handel.mat')
% % whos
%
% audiowrite('Handel.wav',y,Fs)
%
% file_path = dir('Handel.wav');
Uz = unzip('audio_2024-05-12_17-59-56.zip')
file_path = Uz{1}
[audio_data, sample_rate] = audioread(file_path);
SzY = size(audio_data)
Fs = sample_rate
audio_data = audio_data(:,1);
d = audio_data;
% Read an audio waveform
% Calculate the chroma matrix. Use a long FFT to discriminate
% spectral lines as well as possible (2048 is the default value)
cfftlen=2048;
C = chromagram_IF(audio_data,sample_rate,cfftlen);
% The frame advance is always one quarter of the FFT length. Thus,
% the columns of C are at timebase of fftlen/4/sr
tt = (1:size(C,2))*cfftlen/4/sample_rate;
% Plot spectrogram using a shorter window
subplot(311)
sfftlen = 512;
specgram(audio_data,sfftlen,sample_rate);
% Always use a 60 dB colormap range
clim(max(clim)+[-60 0])
% .. and look only at the bottom 4 kHz of spectrum
axis([0 length(d)/sample_rate 0 4000])
title('Original Sound')
% Now the chromagram, also on a dB magnitude scale
subplot(312)
imagesc(tt,1:12,20*log10(C+eps));
axis xy
clim(max(clim)+[-60 0])
title('Chromagram')
function Y = chromagram_IF(d,sr,fftlen,nbin,f_ctr,f_sd)
% Y = chromagram_IF(d,sr,fftlen,nbin,f_ctr,f_sd)
% Calculate a "chromagram" of the sound in d (at sampling rate sr)
% Use windows of fftlen points, hopped by ffthop points
% Divide the octave into nbin steps
% Weight with center frequency f_ctr (in Hz) and gaussian SD f_sd
% (in octaves)
% Use instantaneous frequency to keep only real harmonics.
% 2006-09-26 dpwe@ee.columbia.edu
% Copyright (c) 2006 Columbia University.
%
% This file is part of LabROSA-coversongID
%
% LabROSA-coversongID is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License version 2 as
% published by the Free Software Foundation.
%
% LabROSA-coversongID is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
% General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with LabROSA-coversongID; if not, write to the Free Software
% Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA
% 02110-1301 USA
%
% See the file "COPYING" for the text of the license.
if nargin < 3; fftlen = 2048; end
if nargin < 4; nbin = 12; end
if nargin < 5; f_ctr = 1000; end
if nargin < 6; f_sd = 1; end
A0 = 27.5; % Hz
A440 = 440; % Hz
f_ctr_log = log(f_ctr/A0) / log(2);
fminl = octs2hz(hz2octs(f_ctr)-2*f_sd);
fminu = octs2hz(hz2octs(f_ctr)-f_sd);
fmaxl = octs2hz(hz2octs(f_ctr)+f_sd);
fmaxu = octs2hz(hz2octs(f_ctr)+2*f_sd);
ffthop = fftlen/4;
nchr = 12;
% Calculate spectrogram and IF gram pitch tracks...
[p,m]=ifptrack(d,fftlen,sr,fminl,fminu,fmaxl,fmaxu);
[nbins,ncols] = size(p);
%disp(['ncols = ',num2str(ncols)]);
% chroma-quantized IF sinusoids
Pocts = hz2octs(p+(p==0));
Pocts(p(:)==0) = 0;
% Figure best tuning alignment
nzp = find(p(:)>0);
%hist(nchr*Pmapo(nzp)-round(nchr*Pmapo(nzp)),100)
[hn,hx] = hist(nchr*Pocts(nzp)-round(nchr*Pocts(nzp)),100);
centsoff = hx(find(hn == max(hn)));
% Adjust tunings to align better with chroma
Pocts(nzp) = Pocts(nzp) - centsoff(1)/nchr;
% Quantize to chroma bins
PoctsQ = Pocts;
PoctsQ(nzp) = round(nchr*Pocts(nzp))/nchr;
% map IF pitches to chroma bins
Pmapc = round(nchr*(PoctsQ - floor(PoctsQ)));
Pmapc(p(:) == 0) = -1;
Pmapc(Pmapc(:) == nchr) = 0;
Y = zeros(nchr,ncols);
for t = 1:ncols;
Y(:,t)=(repmat([0:(nchr-1)]',1,size(Pmapc,1))==repmat(Pmapc(:,t)',nchr,1))*m(:,t);
end
end
function hz = octs2hz(octs,A440)
% hz = octs2hz(octs,A440)
% Convert a real-number octave
% into a frequency in Hzfrequency in Hz into a real number counting
% the octaves above A0. So hz2octs(440) = 4.0.
% Optional A440 specifies the Hz to be treated as middle A (default 440).
% 2006-06-29 dpwe@ee.columbia.edu for fft2chromamx
if nargin < 2; A440 = 440; end
% A4 = A440 = 440 Hz, so A0 = 440/16 Hz
hz = (A440/16).*(2.^octs);
end
function octs = hz2octs(freq, A440)
% octs = hz2octs(freq, A440)
% Convert a frequency in Hz into a real number counting
% the octaves above A0. So hz2octs(440) = 4.0
% Optional A440 specifies the Hz to be treated as middle A (default 440).
% 2006-06-29 dpwe@ee.columbia.edu for fft2chromamx
if nargin < 2; A440 = 440; end
% A4 = A440 = 440 Hz, so A0 = 440/16 Hz
octs = log(freq./(A440/16))./log(2);
end
function [p,m,S] = ifptrack(d,w,sr,fminl,fminu,fmaxl,fmaxu)
% [p,m,S] = ifptrack(d,w,sr,fminl,fminu,fmaxl,fmaxu)
% Pitch track based on inst freq.
% Look for adjacent bins with same inst freq.
% d is the input waveform. sr is its sample rate
% w is the basic STFT DFT length (window is half, hop is 1/4)
% S returns the underlying complex STFT.
% fmin,fmax define ramps at edge of sensitivity
% 2006-05-03 dpwe@ee.columbia.edu
% Copyright (c) 2006 Columbia University.
%
% This file is part of LabROSA-coversongID
%
% LabROSA-coversongID is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License version 2 as
% published by the Free Software Foundation.
%
% LabROSA-coversongID is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
% General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with LabROSA-coversongID; if not, write to the Free Software
% Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA
% 02110-1301 USA
%
% See the file "COPYING" for the text of the license.
% downweight fundamentals below here
if nargin < 4; fminl = 150; end
if nargin < 5; fminu = 300; end
% highest frequency we look to
if nargin < 6; fmaxl = 2000; end
if nargin < 7; fmaxu = 4000; end
% Calculate the inst freq gram
[I,S] = ifgram(d,w,w/2,w/4,sr);
% Only look at bins up to 2 kHz
maxbin = round(fmaxu * (w/sr) );
%maxbin = size(I,1)
minbin = round(fminl * (w/sr) );
% Find plateaus in ifgram - stretches where delta IF is < thr
ddif = [I(2:maxbin, :);I(maxbin,:)] - [I(1,:);I(1:(maxbin-1),:)];
% expected increment per bin = sr/w, threshold at 3/4 that
dgood = abs(ddif) < .75*sr/w;
% delete any single bins (both above and below are zero);
dgood = dgood .* ([dgood(2:maxbin,:);dgood(maxbin,:)] > 0 | [dgood(1,:);dgood(1:(maxbin-1),:)] > 0);
% check it out
%p = dgood;
% reconstruct just pitchy cells?
%r = istft(p.*S,w,w/2,w/4);
p = 0*dgood;
m = 0*dgood;
% For each frame, extract all harmonic freqs & magnitudes
for t = 1:size(I,2)
ds = dgood(:,t)';
lds = length(ds);
% find nonzero regions in this vector
st = find(([0,ds(1:(lds-1))]==0) & (ds > 0));
en = find((ds > 0) & ([ds(2:lds),0] == 0));
npks = length(st);
frqs = zeros(1,npks);
mags = zeros(1,npks);
for i = 1:length(st)
bump = abs(S(st(i):en(i),t));
frqs(i) = (bump'*I(st(i):en(i),t))/(sum(bump)+(sum(bump)==0));
mags(i) = sum(bump);
if frqs(i) > fmaxu
mags(i) = 0;
frqs(i) = 0;
elseif frqs(i) > fmaxl
mags(i) = mags(i) * max(0, (fmaxu - frqs(i))/(fmaxu-fmaxl));
end
% downweight magnitudes below? 200 Hz
if frqs(i) < fminl
mags(i) = 0;
frqs(i) = 0;
elseif frqs(i) < fminu
% 1 octave fade-out
mags(i) = mags(i) * (frqs(i) - fminl)/(fminu-fminl);
end
if frqs(i) < 0
mags(i) = 0;
frqs(i) = 0;
end
end
% then just keep the largest at each frame (for now)
% [v,ix] = max(mags);
% p(t) = frqs(ix);
% m(t) = mags(ix);
% No, keep them all
%bin = st;
bin = round((st+en)/2);
p(bin,t) = frqs;
m(bin,t) = mags;
end
%% Pull out the max in each column
%[mm,ix] = max(m);
%% idiom to retrieve different element from each column
%[nr,nc] = size(p);
%pp = p((nr*[0:(nc-1)])+ix);
%mm = m((nr*[0:(nc-1)])+ix);
% r = synthtrax(pp,mm,sr,w/4);
%p = pp;
%m = mm;
end
function [F,D] = ifgram(X, N, W, H, SR)
% [F,D] = ifgram(X, N, W, H, SR) Instantaneous frequency by phase deriv.
% X is a 1-D signal. Process with N-point FFTs applying a W-point
% window, stepping by H points; return (N/2)+1 channels with the
% instantaneous frequency (as a proportion of the sampling rate)
% obtained as the time-derivative of the phase of the complex spectrum
% as described by Toshihiro Abe et al in ICASSP'95, Eurospeech'97
% Same arguments and some common code as dpwebox/stft.m.
% Calculates regular STFT as side effect - returned in D.
% after 1998may02 dpwe@icsi.berkeley.edu
% 2001-03-05 dpwe@ee.columbia.edu revised version
% 2001-12-13 dpwe@ee.columbia.edu Fixed to work when N != W
% $Header: $
% Copyright (c) 2006 Columbia University.
%
% This file is part of LabROSA-coversongID
%
% LabROSA-coversongID is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License version 2 as
% published by the Free Software Foundation.
%
% LabROSA-coversongID is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
% General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with LabROSA-coversongID; if not, write to the Free Software
% Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA
% 02110-1301 USA
%
% See the file "COPYING" for the text of the license.
if nargin < 2; N = 256; end
if nargin < 3; W = N; end
if nargin < 4; H = W/2; end
if nargin < 5; SR = 1; end
s = length(X);
% Make sure it's a single row
if size(X,1) > 1
X = X';
end
%win = [0,hanning(W-1)'];
win = 0.5*(1-cos([0:(W-1)]/W*2*pi));
% Window for discrete differentiation
T = W/SR;
dwin = -pi / T * sin([0:(W-1)]/W*2*pi);
% sum(win) takes out integration due to window, 2 compensates for neg frq
norm = 2/sum(win);
% How many complete windows?
nhops = 1 + floor((s - W)/H);
F = zeros(1 + N/2, nhops);
D = zeros(1 + N/2, nhops);
nmw1 = floor( (N-W)/2 );
nmw2 = N-W - nmw1;
ww = 2*pi*[0:(N-1)]*SR/N;
for h = 1:nhops
u = X((h-1)*H + [1:W]);
% if(h==0)
% plot(u)
% end
% Apply windows now, while the length is right
wu = win.*u;
du = dwin.*u;
% Pad or truncate samples if N != W
if N > W
wu = [zeros(1,nmw1),wu,zeros(1,nmw2)];
du = [zeros(1,nmw1),du,zeros(1,nmw2)];
end
if N < W
wu = wu(-nmw1+[1:N]);
du = du(-nmw1+[1:N]);
end
% FFTs of straight samples plus differential-weighted ones
t1 = fft(fftshift(du));
t2 = fft(fftshift(wu));
% Scale down to factor out length & window effects
D(:,h) = t2(1:(1 + N/2))'*norm;
% Calculate instantaneous frequency from phase of differential spectrum
t = t1 + j*(ww.*t2);
a = real(t2);
b = imag(t2);
da = real(t);
db = imag(t);
instf = (1/(2*pi))*(a.*db - b.*da)./((a.*a + b.*b)+(abs(t2)==0));
% 1/2pi converts rad/s into cycles/s
% sampling rate already factored in as constant in dwin & ww
% so result is in Hz
F(:,h) = instf(1:(1 + N/2))';
end;
end
function X = synthtrax(F, M, SR, SUBF, DUR)
% X = synthtrax(F, M, SR, SUBF, DUR) Reconstruct a sound from track rep'n.
% Each row of F and M contains a series of frequency and magnitude
% samples for a particular track. These will be remodulated and
% overlaid into the output sound X which will run at sample rate SR,
% although the columns in F and M are subsampled from that rate by
% a factor SUBF (default 128). If DUR is nonzero, X will be padded or
% truncated to correspond to just this much time.
% dpwe@icsi.berkeley.edu 1994aug20, 1996aug22
if(nargin<4)
SUBF = 128;
end
if(nargin<5)
DUR = 0;
end
rows = size(F,1);
cols = size(F,2);
opsamps = round(DUR*SR);
if(DUR == 0)
opsamps = 1 + ((cols-1)*SUBF);
end
X = zeros(1, opsamps);
for row = 1:rows
% fprintf(1, 'row %d.. \n', row);
mm = M(row,:);
ff = F(row,:);
% Where mm = 0, ff is undefined. But interp will care, so find points
% and set.
% First, find onsets - points where mm goes from zero (or NaN) to nzero
% Before that, even, set all nan values of mm to zero
mm(find(isnan(mm))) = zeros(1, sum(isnan(mm)));
ff(find(isnan(ff))) = zeros(1, sum(isnan(ff)));
nzv = find(mm);
firstcol = min(nzv);
lastcol = max(nzv);
% for speed, chop off regions of initial and final zero magnitude -
% but want to include one zero from each end if they are there
zz = [max(1, firstcol-1):min(cols,lastcol+1)];
if length(zz) > 0
mm = mm(zz);
ff = ff(zz);
nzcols = prod(size(zz));
mz = (mm==0);
mask = mz & (0==[mz(2:nzcols),1]);
ff = ff.*(1-mask) + mask.*[ff(2:nzcols),0];
% Do offsets too
mask = mz & (0==[1,mz(1:(nzcols-1))]);
ff = ff.*(1-mask) + mask.*[0,ff(1:(nzcols-1))];
% Ok. Can interpolate now
% This is actually the slow part
% % these parameters to interp make it do linear interpolation
% ff = interp(ff, SUBF, 1, 0.001);
% mm = interp(mm, SUBF, 1, 0.001);
% % chop off past-the-end vals from interp
% ff = ff(1:((nzcols-1)*SUBF)+1);
% mm = mm(1:((nzcols-1)*SUBF)+1);
% slinterp does linear interpolation, doesn't extrapolate, 4x faster
ff = slinterp(ff, SUBF);
mm = slinterp(mm, SUBF);
% convert frequency to phase values
pp = cumsum(2*pi*ff/SR);
% run the oscillator and apply the magnitude envelope
xx = mm.*cos(pp);
% add it in to the correct place in the array
base = 1+SUBF*(zz(1)-1);
sizex = prod(size(xx));
ww = (base-1)+[1:sizex];
X(ww) = X(ww) + xx;
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Helper function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Y = slinterp(X,F)
% Y = slinterp(X,F) Simple linear-interpolate X by a factor F
% Y will have ((size(X)-1)*F)+1 points i.e. no extrapolation
% dpwe@icsi.berkeley.edu fast, narrow version for SWS
% Do it by rows
sx = prod(size(X));
% Ravel X to a row
X = X(1:sx);
X1 = [X(2:sx),0];
XX = zeros(F, sx);
for i=0:(F-1)
XX((i+1),:) = ((F-i)/F)*X + (i/F)*X1;
end
% Ravel columns of X for output, discard extrapolation at end
Y = XX(1:((sx-1)*F+1));
end
There is apparently supposed to be a third subplot, however I cannot determine where it is coded.
.
