%design loop filter based on procedure in Karjalainen paper
fs = 44100;
%B = wavread('b.wav');
B = wavread('Bstring 7thfret.wav');
%differentiate B to smooth out signal
B = diff(B);
[b f t] = specgram(B, 2048, fs, blackman(2048), 1024);
figure(1), subplot(311), plot(f,abs(b));
title('STFT of guitar sample');
xlabel('Frequency (Hz)');
ylabel('Magnitude');
axis([0 1e4 0 20]);
%fundamental freq
F=370;
%number of freqs to plot
Nf = 10;
%compute stft:
freqs = F:F:F*Nf;
%freqs = [F/3 F/2 freqs];
window = 1024;
overlap = window/2; %50% overlap
[b f t] = specgram(B, freqs, fs, blackman(window), overlap);
b=b';
%we're only interested in the magnitude spectrum:
b = abs(b);
figure(1), subplot(312), plot(t, b);
title('Guitar sample harmonics vs. time');
ylabel('Magnitude (log scale)');
xlabel('Time (samples)');
[rows cols] = size(b);
%remove first few samples of b - we're not interested in
%components from the initial initial pluck
deltah = b((10:length(b)),(1:cols));
[rows cols] = size(deltah);
for r = 10:rows
for c = 1:cols
deltah(r,c) = 20*log10(deltah(r,c));
end
end
lines = zeros(rows,cols);
slopes = zeros(1,cols);
for c = 1:cols
%fit to straight line:
l = deltah(:,c);
x = t(10:length(t));
m = polyfit(l,x,1);
lines(:,c) = m(1)*x + m(2);
slopes(c) = m(1);
end;
firresp = zeros(1, length(slopes));
for h = 1:length(slopes)
Bk = slopes(h);
%length of delay line
L = F/fs;
%hop size used in stft
H = overlap/fs;
firresp(h) = 10^(Bk*L/(20*H));
end
%assume harmonics after Nfth harmonic decay linearly
step=2*pi*F/fs;
first=2*pi*(F*(Nf+1))/fs;
i=first:step:pi;
%firfill=(.35*(pi-i));
%tweaked firfill to get reasonable attenuation at high frequencies
%(if it decays all the way to zero, the higher harmonics attenuate
%very rapidly)
firfill=.95:-.35/(length(i)-1):.6;
figure(1), subplot(313), stem(freqs,firresp)
title('Extrapolated FIR filter coefficients based on harmonic decay');
xlabel('Frequency (Hz)');
axis([0 freqs(length(freqs)) .9 1]);
[B, A] = invfreqz([firresp firfill], [2*pi*freqs/fs i], 1, 1)
%add a gain factor to keep fundamental freq (low freq) from attenuating
%too rapidly we needed to tweak this filter by hand to get the best results
B= 1.105*B;
figure(2),freqz(B, A);
title('Loop filter response');