How to calculate average frequency and a frequency over time graph all in Matlab?

Hello everyone,

I am encountering an issue with calculating the average frequency in MATLAB. The average frequency result is unexpectedly high and varies with each execution, even though the audio input remains unchanged. Conversely, using the same audio in Python yields more consistent and reasonable results.

If anyone can identify the problem and help fix my MATLAB code, it would be greatly appreciated. Below are two versions of the code: the first records sound directly in MATLAB, and the second processes an existing audio file.

Thank you for your assistance.

Best regards,

UMO

Code 1: also can't figure out freq over time graph

clc
clear
warning off
Fs = 8000; % Sampling frequency in hertz
ch = 1; % Number of channels--2 options-->1 (mono) or 2 (stereo)
datatype = 'uint8';
nbits = 16; % 8, 16, or 24
Nseconds = 5;


% To record audio data from an input device ...
% ... such as a microphone for processing in MATLAB
recorder1 = audiorecorder(Fs, nbits, ch, 1);  % Use device ID 1 for first microphone

No audio input device found on this system.

% Record audio to audiorecorder object,...
% ... hold control until recording completes
recordblocking(recorder1, Nseconds);

disp('End of Recording.');

% Store recorded audio
x1 = getaudiodata(recorder1);

% Write audio file
audiowrite('test1.wav', x1, Fs);

% Calculate time vector
t1 = 0:1/Fs:(length(x1)-1)/Fs;

% Plot time domain signal
figure;
subplot(2, 1, 1);
plot(t1, x1, 'LineWidth', 1.5);
xlabel('Time (sec)');
ylabel('Amplitude');
title('Time Domain Plot of the Recorded Signal');

% Stat Information
maxValue = max(x1);
minValue = min(x1);
meanValue = mean(x1);
stdValue = std(x1);
rmsValue = rms(x1); % root mean square

% Display statistical information
disp('Statistical Information:')
disp(['Max Value: ', num2str(maxValue)])
disp(['Min Value: ', num2str(minValue)])
disp(['Mean Value: ', num2str(meanValue)])
disp(['Standard Deviation: ', num2str(stdValue)])
disp(['RMS Value: ', num2str(rmsValue)])


% Plot time domain signal
%figure;
% subplot(2, 1, 2);
%plot(t2, x2, 'LineWidth', 1.5);
%xlabel('Time (sec)');
%ylabel('Amplitude');
%title('Time Domain Plot of the Recorded Signal mirco 2');

% Calculate FFt
n1 = length(x1);
Y1 = fft(x1);
F1 = Fs/n1*(0:n1-1);

% Find the average frequency

freq = meanfreq( x1 ,Fs );
disp(['Average frequency of the recorded sound: ', num2str(freq), ' Hz']);



% Plot frequency domain signal
% subplot(2, 1, 2);
figure;
plot(F1, abs(Y1(:,1)), 'LineWidth', 1.5);
xlabel('Frequency (Hz)');
ylabel('Amplitude');
title('Frequency Domain Plot of the Audio Signal');

Code 2:

clc
clear
close all
warning off


% Reading input signal
[file, path] = uigetfile('*.wav', 'Select a WAV file');
filename = fullfile(path, file);
[input_signal, Fs] = audioread(filename);

% Create time vector
t = (0:length(input_signal)-1) / Fs;

% Plotting time-domain signal
figure;
plot(t, input_signal);
grid on
title('Time Domain Signal', 'FontName', 'Times New Roman', 'FontSize', 14);
xlabel('Time (sec)', 'FontName', 'Times New Roman', 'FontSize', 12);
ylabel('Amplitude', 'FontName', 'Times New Roman', 'FontSize', 12);

% Stat Information
maxValue = max(input_signal);
minValue = min(input_signal);
meanValue = mean(input_signal);
stdValue = std(input_signal);
rmsValue = rms(input_signal); % root mean square

% Display statistical information
disp('Statistical Information:')
disp(['Max Value: ', num2str(maxValue)])
disp(['Min Value: ', num2str(minValue)])
disp(['Mean Value: ', num2str(meanValue)])
disp(['Standard Deviation: ', num2str(stdValue)])
disp(['RMS Value: ', num2str(rmsValue)])

% Compute FFT
N = length(input_signal);
v = fft(input_signal);
v = v(1:N/2+1); % Keep only the positive frequencies
f = (0:N/2) * Fs / N; % Frequency vector


% Find the average frequency
freq = meanfreq(input_signal, Fs);
disp(['Average frequency of the recorded sound: ', num2str(freq), ' Hz']);

% find avg frequency manually 

tfreq = sum(abs(input_signal).^2); 

% Display the result
disp(['Total frequency domain energy: ', num2str(tfreq)])

[Star Strider]

If you have the Signal Processing Toolbox, use the meanfreq function.

[Ugonne]

Yes I used that install. So I can use meanfreq.

Answer 1

[William Rose]

@Ugonne,

[edit: I adjusted my discussion at the bottom slightly. No changes to code or plots.]

I made a recording with my computer's built-in microphone of the notes C-D-E-F-G-A played on a piano. The piano is reasonably in tune, as you can confirm, if you listen to the recording, and then play the notes on a keyboard.

f=unzip("recordingclip.zip");
[x,fs]=audioread('recordingclip.m4a');
N=length(x); T=N/fs;
fprintf('Sampling rate=%d Hz, duration=%.2f s.\n',fs,T)

Sampling rate=48000 Hz, duration=6.19 s.

This is a stereo recording. Average the L and R channels:

x=mean(x,2);

Find overall mean frequency:

mnFrq=meanfreq(x,fs);
fprintf('Mean frequency=%.1f Hz.\n',mnFrq);

Mean frequency=426.8 Hz.

Find mean frequency in 0.1 s wide window, every 0.05 s (i.e. windows are half-overlapped):

ww=0.1*fs;              % window width (samples)
wOff=floor(ww/2);       % offset of successive windows
nWin=floor(N/wOff-1);   % number of windows
y=zeros(1,nWin);        % allocate array for mean freq
for i=1:nWin
    y(i)=meanfreq(x((i-i)*wOff+1:(i-1)*wOff+ww),fs);
end

Plot mean frequency versus time:

tmf=(1:nWin)*wOff/fs;  % time vector for mean freq (s)
figure; subplot(211)
plot(tmf,y,'-r.');
grid on; xlabel('Time (s)'); ylabel('Mean Freq. (Hz)')

Plot amplitude vs. time:

t=(0:N-1)/fs;
subplot(212); plot(t,x,'-b'); 
grid on; xlabel('Time (s)'); ylabel('Amplitude')

The results are not as expected. Of course it is to be expected that meanfreq() will give strange reuslts when there is hardly any sound, since it is computing (nearly 0) over (nearly 0). But when there is a clearly audible sound, meanfreq() is not always generating the result that my ear hears. The expected frequencies are middle C=262, D=294, E=330, F=349, G=392, A=440. My piano isn't perfect, but it's not that far off. For the first four notes, middle C through F, the plotted mean frequency never approaches the expected value. The plotted mean freq does not change from D to E, even though the change is clear to the listener. For G, and A, the mean frequency plotted above appears correct, eventually, although I am surprised at how long it take to settle - it seems to take longer than the window width.

The difference between what meanfreq() is reporting and what my ear hears is due to a variety of neurophysiological factors, having to do with how the ear weights intensities (logarithmically); the fact that we tend to identify the mode of the spectrum, not the mean, when we recognize pitch; nonlinearities; etc.

[William Rose]

@Ugonne,

I notice that your code includes a section at the end where you start to try to compute mean frequency "manually". Here is a way to find the mean frequency without using meanfreq(). We will make two signals to demonstrate.

fs=8000;
N=5*fs;                 % number of samples for 5 seconds of data
t=(0:N-1)/fs;           % time (s)
x1=cos(2*pi*1000*t);    % x1=pure tone at 1000 Hz
x2=randn(1,N);          % x2=white noise

Find the mean frequency of each signal:

f=(0:N/2)*fs/N;                 % frequencies 0 through fs/2 (Hz)
Pxx1=abs(fft(x1)).^2;
Pxx1=Pxx1(1:N/2+1);             % power 0 through fs/2
Pxx2=abs(fft(x2)).^2;
Pxx2=Pxx2(1:N/2+1);             % power 0 through fs/2
mnFrq1=sum(Pxx1.*f)/sum(Pxx1);  % mean freq (Hz)
mnFrq2=sum(Pxx2.*f)/sum(Pxx2);  % mean freq (Hz)

Display the mean frequencies found manually:

fprintf('x1, x2 mean frequency (Hz)=%.1f, %.1f.\n',mnFrq1,mnFrq2)

x1, x2 mean frequency (Hz)=1000.0, 2004.3.

x1 is a pure tone at 1000 Hz, so the result for x1 is as expected.

x2 is white noise from 0 to fs/2=4 kHz, so the expected mean frequency is 2 kHz. The actual value varies somewhat due to the random nature of x2.