How to plot Chroma_Features for an audio file?
5 views (last 30 days)
Show older comments
%% Read in the file
clearvars;
close all;
file_path = 'C:/Users/nimae/Downloads/audio_2024-05-12_17-59-56.ogg';
[audio_data, sample_rate] = audioread(file_path);
%% Play original file
% pOrig = audioplayer(audio_data, sample_rate);
% pOrig.play;
%% Plot both audio channels
N = size(audio_data,1); % Determine total number of samples in audio file
figure;
subplot(2,1,1);
stem(1:N, audio_data(:,1));
title('Left Channel');
subplot(2,1,2);
% Increase resolution by increasing the number of points in FFT
N = 2^nextpow2(length(audio_data(:,1)));
% Plot the spectrum with increased resolution
df = sample_rate / N;
w = (-N/2 : N/2 - 1) * df;
y = fft(audio_data(:,1), N) / N; % For normalizing, but not needed for our analysis
y2 = fftshift(y);
% Define the desired frequency range
f_min = 0; % Minimum frequency in Hz
f_max = 5.6e4; % Maximum frequency in Hz
% Find the corresponding indices in the frequency axis
idx_min = find(w >= f_min, 1, 'first');
idx_max = find(w <= f_max, 1, 'last');
% Plot the spectrum with increased resolution and within the desired frequency range
figure;
plot(w(idx_min:idx_max), abs(y2(idx_min:idx_max)));
xlabel('Frequency (Hz)');
ylabel('Magnitude');
title('Spectrum of Audio Signal');
grid on;
hold on; % Add subsequent plots to the same figure
% Frequencies you provided
frequencies = [
16.35 32.702 65.404 130.808 261.616 523.232 1046.464 2092.928 4185.856 8371.712;
17.32 34.648 69.296 138.592 277.184 554.368 1108.736 2217.472 4434.944 8869.888;
18.35 36.708 73.416 146.832 293.664 587.328 1174.656 2349.312 4698.624 9397.248;
19.45 38.89 77.78 155.56 311.12 622.24 1244.48 2488.96 4977.92 9955.84;
20.60 41.202 82.404 164.808 329.616 659.232 1318.464 2636.928 5273.856 10547.712;
21.83 43.654 87.308 174.616 349.232 698.464 1396.928 2793.856 5587.712 11175.424;
23.12 46.248 92.496 184.992 369.984 739.968 1479.936 2959.872 5919.744 11839.488;
24.50 48.998 97.996 195.992 391.984 783.968 1567.936 3135.872 6271.744 12543.488;
25.96 51.912 103.824 207.648 415.296 830.592 1661.184 3322.368 6644.736 13289.472;
27.50 55 110 220 440 880 1760 3520 7040 14080;
29.14 58.27 116.54 233.08 466.16 932.32 1864.64 3729.28 7458.56 14917.12;
30.87 61.736 123.472 246.944 493.888 987.776 1975.552 3951.104 7902.208 15804.416
];
% Plot the frequencies you provided
for i = 1:size(frequencies, 1)
plot(frequencies(i, :), zeros(size(frequencies, 2)), 'o', 'MarkerSize', 5);
end
legend('Spectrum', 'Frequencies provided');
% Convert frequencies to chroma values
chroma_values = mod(round(log2(frequencies/440) * 12), 12);
% Plot the chroma diagram
figure;
imagesc(chroma_values);
colormap(jet);
colorbar;
xlabel('Time');
ylabel('Chroma');
title('Chroma Diagram');
function Y = chromagram_IF(d,sample_rate,fftlen,~,f_ctr,f_sd)
%% Function Definitions
% 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')
% 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);
[~,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] = histogram(nchr*Pocts(nzp)-round(nchr*Pocts(nzp)),100);
centsoff = hx(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
hello,
i have coded this file and copied most of it from Columbian University, as you can see. I want to plot simply the chroma features for my file in variation of time, i mean determining exactly at which time which musical note is played in the file. For that we need at the X axis, time, and in Yaxis Frequencies of those musical notes(which i have provided in the "Frequencies" Matrix).
So I'd be thankful if someone can guide me a little bit through this. My code runs halfway but i don't get the desired Chroma Diagram, said above.
0 Comments
Accepted Answer
Star Strider
on 13 May 2024
Edited: Star Strider
on 13 May 2024
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.
.
8 Comments
Star Strider
on 14 May 2024
As always, my pleasure!
Thank you! Yours is an interesting question, and I learned an aspect of signal processing that I never before knew existed by answering it.
More Answers (0)
See Also
Categories
Find more on Acoustics, Noise and Vibration in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!