Finding MNF for every second inteval

Hello everyone! I'm currently trying to plot the MNF value from every 1 second interval of filtered signal. Here is the code I write:

%% MNF and MDF of signal
%Calculate MNF from every 1 s interval of the filtered signal
start = 0;
stop = 120;
%num_windows = floor(length(psd1) / fs);
interval=ceil((stop-start)/fs);
t_mnf = linspace(start,stop,interval+1);

mnf_s = zeros(1, interval+1);

for i=start:stop
    % Check for starting index (avoid going below 1)
    start_index = max(1, (i-1)*fs+1);
    % Check for ending index (avoid exceeding psd1 length)
    end_index = min(length(psd1), i*fs);
    %Extract current window data by assuming 1 second window
    window_data = psd1(start_index:end_index);
    mnf_s(i-start+1)=meanfreq(window_data,fs);
end

figure(5)
plot(t,mnf_s,'-*k')
xlabel('Time (s)'); ylabel('MNF (Hz)')
title('MNF')

The 'psd1' variable contains the PSD value of filtered signal. The code doesn't work. Where could it went wrong?

[Mathieu NOE]

hello

it would help if you could supply the data as well

but I have a feeling the problem is here

mnf_s(i-start+1)=meanfreq(window_data,fs);

^^

according to meanfreq doc

freq = meanfreq(x,fs) estimates the mean frequency in terms of the sample rate, fs.

freq = meanfreq(pxx,f) returns the mean frequency of a power spectral density (PSD) estimate, pxx. The frequencies, f, correspond to the estimates in pxx.

what you have coded is kinda mix up of both possibilities , but that is not gonna work

please provide the frequency vector (and not the sampling frequency) when dealing with psd input

[Mathieu NOE]

BTW, the code does not show how you generate the 1s data PSD - did you remove that from the posted code intentionnaly ?

[Keisha Alfreda]

Thank you so much for the respond. And this is how I generate the PSD data

%% Find Power Spectral Density from FFT
freq = 0:fs/L:fs-1/fs;

vv_psd = abs(fft(vv_filter)).^2/L; %Two-sided power from filtered signal
psd1 = vv_psd(1:L/2+1); %One-sided power
psd1(2:end-1) = 2*psd1(2:end-1);

figure(4)
plot(freq(1:L/2+1),psd1)
xlabel('Frequency (Hz)')
xlim([0 fs/2])
ylabel('Power density (dB)/Hz')
title('Power Spectral Density using FFT')

This is the result of the one-sided PSD

[Mathieu NOE]

NB your ylabel says psd is in dB / Hz but you are plotting a psd in linear units ² / Hz

[Mathieu NOE]

it's not clear for me , but your main code should compute the psd for every second of data

is it what you are doing ?

Answer 1

[Sudarsanan A K]

Hi Keisha,

The provided code snippet aims to calculate the Mean Frequency (MNF) from every 1-second interval of a filtered signal represented by its Power Spectral Density (PSD) data, "psd1". However, there are a few issues and potential improvements in the implementation:

• Incorrect Loop Range: The loop iterates from "start" to 'stop" (0 to 120 in this case), which is likely meant to represent seconds. However, the loop variable "i" is used as if it represents indices directly. Since "fs" (sampling frequency) is involved in calculating indices, the loop should iterate through intervals, not through every integer from "start" to "stop".
• Misinterpretation of "fs" with "meanfreq": The "meanfreq" function expects a vector of frequencies corresponding to the PSD data points, not the sampling frequency. The second argument should be a frequency vector that matches the PSD data in "window_data".

Also, you could use "pwelch" function to easily compute the PSD and the frequency vector.

Here is a simple demonstration how you could compute the MNF:

% Define the sampling frequency and time vector
fs = 1000; % Sampling frequency in Hz
maxTime = 5; % In seconds
t = 0:1/fs:maxTime-1/fs; % 5 seconds duration

% Generate a linear chirp signal (Sample signal)
f0 = 1; % Start frequency of the chirp in Hz
f1 = 50; % End frequency of the chirp in Hz
signal = chirp(t, f0, t(end), f1);

% Calculate the Mean Frequency (MNF) for every 1-second interval of the signal
interval = 1; % 1 second interval
num_windows = floor(length(t) / fs / interval);
mnf_s = zeros(1, num_windows);
t_mnf = linspace(0, num_windows-1, num_windows);

for i = 1:num_windows
    start_index = (i-1) * fs + 1;
    end_index = i * fs;
    % Ensure the end_index does not exceed the length of the signal
    end_index = min(length(signal), end_index);
    % Extract current window data
    window_data = signal(start_index:end_index);
    % Calculate MNF for the current window using the pwelch method to estimate PSD
    [Pxx, F] = pwelch(window_data, 500, 300, 500, fs);
    mnf_s(i) = meanfreq(Pxx, F);
end

% Plot the chirp signal
figure(1);
plot(t, signal);
title('Chirp Signal');
xlabel('Time (s)');
ylabel('Amplitude');

% Plot the Mean Frequency over time
figure(2);
plot(t_mnf, mnf_s, '-*k');
xlabel('Time (s)');
ylabel('MNF (Hz)');
title('Mean Frequency (MNF) of the Chirp Signal');

% Display the calculated MNF values
disp('Calculated MNF values for each 1-second interval:');

Calculated MNF values for each 1-second interval:

disp(mnf_s);

    5.4394   15.2357   25.0337   34.8293   44.6166

Explanation:

• Chirp Signal Generation: The "chirp" function generates a signal whose frequency increases linearly from "f0" Hz to "f1" Hz over "maxTime" seconds. This signal is a perfect test case for observing how the MNF calculation adapts to changing frequencies.
• MNF Calculation: The code calculates the MNF for each 1-second interval of the chirp signal. Due to the linear frequency increase of the chirp, we expect the MNF to show a clear trend that reflects this change over each interval.
• Plotting and Output: The code plots both the chirp signal and the MNF values over time. The plot of MNF values should show an increase, reflecting the frequency increase in the chirp signal.

The following resources might be useful for further understanding of the calculations of PSD and MNF using MATLAB:

• "meanfreq": https://www.mathworks.com/help/signal/ref/meanfreq.html
• "pwelch": https://www.mathworks.com/help/signal/ref/pwelch.html

I hope this helps!

[Keisha Alfreda]

Thank you very much for the answer, helps me a lot to learn better!