Can someone tell me how I can implement Harmonic Product Spectrum using MATLAB to find the fundamental frequency of a note in the presence of harmonics?? I know I'm supposed to downsample my signal a number of times (after performing fft of course) and then multiply them with the original signal.
Say my fft signal is "FFT1"
then the code would roughly be like
hps1 = downsample(FFT1,2);
hps2 = downsample(FFT1,3);
hps = FFT1.*hps1.*hps2;
Is this code correct??? I want to know if I've downsampled properly and since each variable has a different length multiplying them results in matrix dimension error.. I really need some real quick help as its for a project work... Really desperate.... Thanx in advance....
OK you can't do "hps = FFT1.*hps1.*hps2;" for each downsampled data, do you have different sizes ...
I did a example for you how make a very simple Harmonic Product Spectrum (HPS) using 5 harmonics decimation (downsample), I just test in sinusoidal signals, I get very near fundamental frequency in my tests.
This code only shows how to compute the main steps of the algorithm, is very likely that you will need improve it !
Source:
%[x,fs] = wavread('ederwander_IN_250Hz.wav');
CorrectFactor = 0.986;
threshold = 0.2;
%F0 start test
f = 250;
fs = 44100;
signal= 0.9*sin(2*pi*f/fs*(0:9999));
x=signal';
framed = x(1:4096);
windowed = framed .* hann(length(framed));
FFT = fft(windowed, 4096);
FFT = FFT(1 : size(FFT,1) / 2);
FFT = abs(FFT);
hps1 = downsample(FFT,1);
hps2 = downsample(FFT,2);
hps3 = downsample(FFT,3);
hps4 = downsample(FFT,4);
hps5 = downsample(FFT,5);
y = [];
for i=1:length(hps5)
Product = hps1(i) * hps2(i) * hps3(i) * hps4(i) * hps5(i);
y(i) = [Product];
end
[m,n]=findpeaks(y, 'SORTSTR', 'descend');
Maximum = n(1);
%try fix octave error
if (y(n(1)) * 0.5) > (y(n(2))) %& ( ( m(2) / m(1) ) > threshold )
Maximum = n(length(n));
end
F0 = ( (Maximum / 4096) * fs ) * CorrectFactor
plot(y)
HPS usually generates an error showing the pitch one octave up, I change a bit a code, see above :-)