MATLAB script to determine the frequency that has the greatest power

Question

I am doing the linear motion measurement with a moving audio source and a stationary observer. Described here: http://www.animations.physics.unsw.edu.au/labs/Doppler/doppler_lab.html

How to write a MATLAB script to do get the sample number of the frequency that has the greatest power in the audio file?

Answer 1

What you need is the Time-Frequency localization information. This can be obtained using Short-time Fourier transform. There are many other Time-Frequency analysis techniques, STFT being the simplest and hence a good starting point. Here is a simple code to help you understand the concept:

% Parameters
Fs = 44100; % Sample rate (44100 Hz)
t = 0:1/Fs:5; % Time instances
t1 = 5; % End time of signal, 5 secs
f0 = 440; % frequency swiped from 440 Hz
f1 = 880; % to 880 Hz

% Signal generation
audio = chirp(t,f0,t1,f1);
% wavplay(audio,Fs) % to play the audio


% Signal analysis
window = 2050; % Should be minimum twice the maximum frequency we want to analyze
noverlap = 1025; % 50% overlap
nfft = 44100;

[S,F,T,P] = spectrogram(audio,window,noverlap,nfft,Fs); % Spectrogram takes the STFT of the signal
% P matrix contains the power spectral density of each segment of the STFT

% Plotting results
imagesc(real(S)) % frequency-time Plot of the signal
ylim([0 1000])
xlabel('Time (secs)')
ylabel('Frequency (Hz)')
title('Time-Frequency plot of a Audio signal')

To get the sample number, you just need to find the time instance at which the frequency of your interest appears, and use sampling frequency to compute the sample number.

P is the power spectral density matrix. Along y-axis is the frequency, x-axis is time, and power contributed by each frequency at every instant is stored in this matrix. You need the element which has the highest value in the entire matrix. A code like below should work:

[maxElement, maxElementTimeIndex] = max(max(P, [], 1)); % max(P, [], 1) finds maximum power for each time instance, max(max(P,[],1)) finds maximum for the entire 2D matrix.
maxPoweredSampleInAudioSignal = (maxElementTimeIndex-1) * Fs; % This calculation is made within the limitations of STFT, so approximately here the maximum power for any frequency is present