1

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....

ederwander
  • 3,410
  • 1
  • 18
  • 23
user2482542
  • 361
  • 2
  • 14
  • 26
  • It's unclear what you're asking. You should be more specific – Luis Mendo Nov 04 '13 at 10:11
  • Well I want to do a Harmonic Product Spectrum analysis on my audio signal after performing STFT so i can get the correct fundamental frequency when there are harmonics present. I'm not familiar with HPS so i want some help in writing a MATLAB code that would implement it... – user2482542 Nov 04 '13 at 10:26

1 Answers1

2

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 :-)

ederwander
  • 3,410
  • 1
  • 18
  • 23
  • Thank you VERY much.. I do have a doubt though. I thought the initial FFT data is downsampled by 2,3 etc.. However in your code you use `hps3 = downsample(hps2,3); ` where you are downsampling the already downsampled signal.. A bit confused there... And why are you using a "correctFactor"? – user2482542 Nov 04 '13 at 16:14
  • take a look here [link](http://cnx.org/content/m11714/latest/HPS.jpg), you can see in this picture abs(FFT) downsampled by /1, /2, /3, /N, YEAH I just tested downsampling at the already downsampled, but you can do it the FFT, the correctFactor is just to increase the results from sinusoidal tests – ederwander Nov 04 '13 at 17:40
  • Ok thank you.. I'll try it both ways and see.. By the way, how exactly did you arrive at the correctfactor value?? Is there a mathematical procedure to follow?? Would the algorithm still work if I didnt use a correctFactor?? – user2482542 Nov 05 '13 at 05:23
  • Hi change a bit the code see above, the correctfactor I only measure by observation, yeah this work without correctfactor, how said for you, I only tested in sinusoids, I believe that you will need improve some things, I did it yesterday just following some papers :-) – ederwander Nov 05 '13 at 14:06
  • Thank you. But sadly enough, even with the corrections you've made, it doesn't seem to work for me :( – user2482542 Nov 07 '13 at 13:28