Sunday, 15 September 2013

Finding prominent frequencies from FFT

Finding prominent frequencies from FFT

Im trying to find the prominent peaks in an audio signal (piano recording)
using STFT. This is what Ive done so far 1. Obtain the envelope of the
time domain signal 2. Determine the peaks in the enveloped signal and use
them as note onsets 3. perform FFT for the samples between each 2
consecutive onsets.
Now that I have the FFT, I want to find the peaks corresponding to the
notes played... when i try using the findpeaks function at some points it
says its an empty matrix. Im not sure how exactly to do this so need some
help urgently...... Thank you
clear all;
clear max;
clc;
[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);
%windowing = hamming(32768); %Windowing function
% Plot time domain signal
figure(1);
subplot(2,1,1)
plot(t,3*song)
title('Wave File')
ylabel('Amplitude')
xlabel('Length (in seconds)')
%ylim([-1.1 1.1])
xlim([0 N/FS])
%----------------------Finding the envelope of the signal-----------------%
% Gaussian Filter
x = 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 = -x .* exp( -(x.^2)/(2*sigma.^2)); % compute first derivative,
but leave constants out
myFilter = myFilter / sum( abs( myFilter ) ); % normalize
% Plot Gaussian Filter
subplot(2,1,2)
plot(myFilter)
title('Edge Detection Filter')
% 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.01
maxtab = [maxtab ; max_pos max1];
min1 = peak;
min_pos = x(i);
lookformax = 0;
end
else
if peak > min1+0.05
mintab = [mintab ; min_pos min1];
max1 = peak;
max_pos = x(i);
lookformax = 1;
end
end
end
% % Plot song filtered with edge detector
figure(2)
plot(1/FS:1/FS:N/FS,tedges)
title('Song Filtered With Edge Detector 1')
xlabel('Time (s)')
ylabel('Amplitude')
ylim([-1 1.1])
xlim([0 N/FS])
hold on;
plot(maxtab(:,1)/FS, maxtab(:,2), 'ro')
plot(mintab(:,1)/FS, mintab(:,2), 'ko')
max_col = maxtab(:,1);
peaks_det = max_col/FS;
No_of_peaks = length(peaks_det);
song = detrend(song);
%---------------------------Performing FFT--------------------------------%
for i = 2:No_of_peaks
song_seg = song(max_col(i-1):max_col(i)-1);
% song_seg = song(max_col(6):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;
N=5;Fst1=50;Fp1=60; Fp2=1040; Fst2=1050;
% d = fdesign.bandpass('N,Fst1,Fp1,Fp2,Fst2');
% h = design(d);
% seg_fft = filter(h, seg_fft);
% seg_fft(1) = 0;
%
f = FS/2*linspace(0,1,NFFT/2+1);
seg_fft2 = 2*abs(seg_fft(1:NFFT/2+1));
L5 = length(song_seg);
figure(1+i)
plot(f,seg_fft2)
title('Frequency spectrum of signal')
xlabel('Frequency (Hz)')
%xlim([0 2500])
ylabel('|Y(f)|')
ylim([0 300])
%[B, IX] = sort(seg_fft2)
%[points loc] = findpeaks(seg_fft);
%STFT_out(:,i) = seg_fft2;
%P=max(seg_fft2)
[points, loc] = findpeaks(seg_fft2,'THRESHOLD',20)
end

No comments:

Post a Comment