Plotting Harmonic Product Spectrum
10 views (last 30 days)
Show older comments
I used Harmonic Product Spectrum to find the fundamental note present when a number of harmonics are present. This is the code I've implemented;
[song,FS] = wavread('C major.wav');
%sound(song,FS);
P = 20000;
N=length(song); % length of song
t=0:1/FS:(N-1)/FS; % define time period
song = sum(song,2);
song=abs(song);
%----------------------Finding the envelope of the signal-----------------%
% Gaussian Filter
w = linspace( -1, 1, P); % create a vector of P values between -1 and 1 inclusive
sigma = 0.335; % standard deviation used in Gaussian formula
myFilter = -w .* exp( -(w.^2)/(2*sigma.^2)); % compute first derivative, but leave constants out
myFilter = myFilter / sum( abs( myFilter ) ); % normalize
% fft convolution
myFilter = myFilter(:); % create a column vector
song(length(song)+length(myFilter)-1) = 0; %zero pad song
myFilter(length(song)) = 0; %zero pad myFilter
edges =ifft(fft(song).*fft(myFilter));
tedges=edges(P:N+P-1); % shift by P/2 so peaks line up w/ edges
tedges=tedges/max(abs(tedges)); % normalize
%---------------------------Onset Detection-------------------------------%
% Finding peaks
maxtab = [];
mintab = [];
x = (1:length(tedges));
min1 = Inf;
max1 = -Inf;
min_pos = NaN;
max_pos = NaN;
lookformax = 1;
for i=1:length(tedges)
peak = tedges(i:i);
if peak > max1,
max1 = peak;
max_pos = x(i);
end
if peak < min1,
min1 = peak;
min_pos = x(i);
end
if lookformax
if peak < max1-0.07
maxtab = [maxtab ; max_pos max1];
min1 = peak;
min_pos = x(i);
lookformax = 0;
end
else
if peak > min1+0.08
mintab = [mintab ; min_pos min1];
max1 = peak;
max_pos = x(i);
lookformax = 1;
end
end
end
max_col = maxtab(:,1);
peaks_det = max_col/FS;
No_of_peaks = length(peaks_det);
[song,FS] = wavread('C major.wav');
song = sum(song,2);
%---------------------------Performing STFT--------------------------------%
h = 1;
%for i = 2:No_of_peaks
song_seg = song(max_col(7-1):max_col(7)-1);
L = length(song_seg);
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
seg_fft = fft(song_seg,NFFT);%/L;
f = FS/2*linspace(0,1,NFFT/2+1);
seg_fft_2 = 2*abs(seg_fft(1:NFFT/2+1));
L5 = length(song_seg);
figure(6)
plot(f,seg_fft_2)
%plot(1:L/2,seg_fft(1:L/2))
title('Frequency spectrum of signal (seg_fft)')
xlabel('Frequency (Hz)')
xlim([0 2500])
ylabel('|Y(f)|')
ylim([0 500])
%----------------Performing Harmonic Product Spectrum---------------------%
% In harmonic prodcut spectrum, you downsample the fft data several times and multiply all those with the original fft data to get the maximum peak.
%HPS
seg_fft = seg_fft(1 : size(seg_fft,1)/2 );
seg_fft = abs(seg_fft);
a = length(seg_fft);
seg_fft2 = ones(size(seg_fft));
seg_fft3 = ones(size(seg_fft));
seg_fft4 = ones(size(seg_fft));
seg_fft5 = ones(size(seg_fft));
for i = 1:((length(seg_fft)-1)/2)
seg_fft2(i,1) = seg_fft(2*i,1);%(seg_fft(2*i,1) + seg_fft((2*i)+1,1))/2;
end
%b= size(seg_fft2)
L1 = length(seg_fft2);
NFFT1 = 2^nextpow2(L1); % Next power of 2 from length of y
f1 = FS/2*linspace(0,1,NFFT1/2+1);
seg_fft12 = 2*abs(seg_fft2(1:NFFT1/2+1));
figure(7);
plot(f1,seg_fft12)
title('Frequency spectrum of signal (seg_fft2)')
xlabel('Frequency (Hz)')
xlim([0 2500])
ylabel('|Y(f)|')
ylim([0 500])
This is the plot for Figure 6
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/153099/image.jpeg)
So in the real scenario once I perform HPS (downsample by 2) the peak at 440.1 should shift down to 220 while the peak at 881 should shift down to around 440. But when I plot the graph that is not what I get. Insetad this is the graph I get,
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/153100/image.jpeg)
Why don't I get the correct graph???? I don't seem to understand what Im doing wrong here... Can someone please have a look and let me know.. Thank you.....
0 Comments
Answers (0)
See Also
Categories
Find more on Time-Frequency Analysis in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!