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?
5 个评论
Mathieu NOE
2024-4-22
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(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
2024-4-22
BTW, the code does not show how you generate the 1s data PSD - did you remove that from the posted code intentionnaly ?
Keisha Alfreda
2024-4-23
编辑：Keisha Alfreda
2024-4-23
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
2024-4-23
NB your ylabel says psd is in dB / Hz but you are plotting a psd in linear units ² / Hz
Mathieu NOE
2024-4-24
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 ?
采纳的回答
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:');
disp(mnf_s);
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!
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Parametric Spectral Estimation 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
