Mean and Median Frequency

Hi, I am a newbie in doing coding for mean and median frequency. first of all, i cut the signal before applying the fft. then i apply window hamming to fft then i create coding for mean and median frequency. i want to extract the value for mean and median but it gives output columns zero. unfortunately i couldnt detect the problem regarding my coding. can someone help me please?

b=EMG_Filtered2segment;               %signal data
fs =1000;                              %sampling frequency
cut=1*fs;                                %cut signal every second
Window=hamming(cut);
L= 300
id1=1
for i = 1:length(b)/cut
    slice=b(1+(i-1)*cut: i*cut);
    signal(:,i)=slice;
    Ls=length(slice);
    Nfftslice=2^nextpow2(Ls);

    signalwin=Window.*slice;
    FFT=fft(signalwin);
    absFFT=abs(FFT);
    FFTonesided=2*absFFT(1:Ls/2+1); %double amplitude, one side
    FFTonesided(1)=FFTonesided(1)/2; %DC part is not doubled
    Hspec=spectrum.welch('Hamming',500,25);
    hspd=psd(Hspec,slice,'NFFT',Nfftslice,'Fs',fs);
    Pw = hpsd.Data;
    Fw=FFTonesided;
    PW=cumsum(Fw);
    Pwdiv2=median(PW);
    Fs=1000;
    Freq= Fs*(0:(Ls/2))/Ls;

end

%MEAN AND MEDIAN%
% a=find(FFTonesided(:,1));
for t=1:length(b)/cut
    meanPSD(:,i) = sum(Freq.*FFTonesided) / sum(FFTonesided);
    OverHalfIdx(:,i) = find(PW>=Pwdiv2,1,'first');
    UnderHalfIdx(:,i)= find(PW<=Pwdiv2,1,'last');
    MedianPSD(:,i) = (Freq(OverHalfIdx(:,i))+Freq(UnderHalfIdx(:,i)))/2
end

[KALYAN ACHARJYA]

EMG_Filtered2segment ??

[ZackyAs]

Hello sir. Hope u doing ok there. Basically, EMG_filtered2segment is my data subject. I attached below the full coding so that u can picture what I want to ask

close all
clear all

L=importdata('BRYAN 1ST HALF (1).txt');
T1_10rm_BS_2=importdata('BRYAN 2ND HALF.txt');
disp(L)
L1=L(:,1)
whos 
n=numel(L);
segment=L(round((n)/30):(n-1/2));
segment=segment(round((n)/3:(2*n)/3))
segment1=T1_10rm_BS_2(round(1:(n-3)));
segment1=segment1(round((n)/3:(2*n)/3))
t1=length(segment)
t11=length(segment1)
N=length(L1)
fs=1000;
Nyquist = fs/2;
t=1/fs                                                                   
T=(0:t1-1)*t;
T1=(0:t11-1)*t;
figure  
plot(T,segment,T1,segment1,'-r');
grid on;
title('raw signals')



%zpk filter for segment 1
[z,p,k] = butter(6,20/(fs/2),'high');
[sos,g] =zp2sos(z,p,k);
EMG_Filteredsegment = filtfilt(sos,g,segment);                  % HP Filter EMG
[z,p,k] = butter(4,450/(fs/2),'low');
[sos,g] =zp2sos(z,p,k);
EMG_Filtered2segment = filtfilt(sos,g,EMG_Filteredsegment);            % LP Filter EMG

figure
grid
subplot (2,1,1)
plot(EMG_Filteredsegment)
ylim([-0.5 0.3])
title ('Partial Filtered Segment 1: Time-Domain')
subplot (2,1,2)
plot(EMG_Filtered2segment)
ylim([-0.5 0.3])
title ('Filtered Signals Segment 1: Time-Domain')


%zpk filter for segment 2
[z,p,k] = butter(6,20/(fs/2),'high');
[sos,g] =zp2sos(z,p,k);
EMG_Filteredsegment1 = filtfilt(sos,g,segment1);                  % HP Filter EMG
[z,p,k] = butter(4,450/(fs/2),'low');
[sos,g] =zp2sos(z,p,k);
EMG_Filtered2segment1 = filtfilt(sos,g,EMG_Filteredsegment1);            % LP Filter EMG

figure(3)
grid
subplot (2,1,1)
plot(EMG_Filteredsegment1)
ylim([-0.5 0.4])
title ('Partial Filtered Segment 2: Time-Domain')
subplot (2,1,2)
plot(EMG_Filtered2segment1)
ylim([-0.5 0.4])
title ('Filtered Signals Segment 2: Time-Domain')

figure
grid
plot(T,EMG_Filtered2segment,T1,EMG_Filteredsegment1,'-r');
ylim([-0.5 0.3])
title ('Filtered Signals Both Segment')

%%RMS 1ST HALF DATA%%

   segments=1;
TS=300; % this is not sample time, Ts is always time difference bewteen 2 measurement points
len = round(segments/TS);
n=length(segment)
L=300
segV=zeros(L,1);
for id =len+1:L
    segV(id) = rms(EMG_Filtered2segment(id-len : id) );
    segV1(id) = rms(EMG_Filtered2segment1(id-len : id) );
    II1=segV(id)
end

%%RMS 2ND HALF DATA%%

segment2=1;
TSS=300
len1 = round(segment2/TSS);
n1=length(segment1)
L=300
segV=zeros(L,1);
for id1 =len1+1:L
    segV1(id1) = rms(EMG_Filtered2segment1(id1-len1 : id1) );
    segV(id1) = rms(EMG_Filtered2segment(id1-len1 : id1) );
end


figure
subplot(2,1,1)
plot(segV)
title('rms of 1st half')
subplot(2,1,2)
plot(segV1, 'r')
title('rms OF 2ND HALF')

b=EMG_Filtered2segment;               %signal data
fs =1000;                              %sampling frequency
cut=1*fs;                                %cut signal every second
Window=hamming(cut);
L= 300
id1=1
for i = 1:length(b)/cut
    slice=b(1+(i-1)*cut: i*cut);
    signal(:,i)=slice;
    Ls=length(slice);
    Nfftslice=2^nextpow2(Ls);

    signalwin=Window.*slice;
    FFT=fft(signalwin);
    absFFT=abs(FFT);
    FFTonesided=2*absFFT(1:Ls/2+1); %double amplitude, one side
    FFTonesided(1)=FFTonesided(1)/2; %DC part is not doubled
    Fw=FFTonesided;
    PW=cumsum(Fw);
    Pwdiv2=median(PW);
    Fs=1000;
    Freq= Fs*(0:(Ls/2))/Ls;

end

%MEAN AND MEDIAN%
% a=find(FFTonesided(:,1));
for t=1:length(b)/cut
    meanPSD(:,i) = sum(Freq.*FFTonesided) / sum(FFTonesided);
    OverHalfIdx(:,i) = find(PW>=Pwdiv2,1,'first');
    UnderHalfIdx(:,i)= find(PW<=Pwdiv2,1,'last');
    MedianPSD(:,i) = (Freq(OverHalfIdx(:,i))+Freq(UnderHalfIdx(:,i)))/2
end

[Image Analyst]

How are we supposed to run these lines when you're not supplying us wtih the files?

L=importdata('BRYAN 1ST HALF (1).txt');
T1_10rm_BS_2=importdata('BRYAN 2ND HALF.txt');

Please read: TUTORIAL: How to ask a question (on Answers) and get a fast answer

[ZackyAs]

I am sorry sir. here

[Image Analyst]

So you have a couple of signals, and there are a bunch of frequencies, and for each frequency there is a relative power at that frequency.

What do you mean by the mean and median frequency? Do you just want the half way frequency between 0 and the maximum frequency? Or do you want the mean weighted by the power at each frequency?

[ZackyAs]

Hello there. Basically, the one i am analysing is changes of muscle fatigue characteristics during driving using emg signal. so my parameters will be time domain which is RMS and also frequency domain which is Mean and Median frequency. so for emg fatigue condition, RMS will increase and Mean and median frequency decrease. to be honest, im still lacking the knowledge for mean and median frequency thats why i ask for what is missing there. can u help me and suggest me sir? i am apologise as i couldnt answer your questions sir since i have no idea about it.

Answer 1

[Star Strider]

If you have R2015a or later, use the medfreq and meanfreq functions.

[ZackyAs]

Hello sir. you meant instead of using this meanpsd n medianpsd, use med n meanfreq functions? it change the whole coding isnt it?

    meanPSD(:,i) = sum(Freq.*FFTonesided) / sum(FFTonesided);
    OverHalfIdx(:,i) = find(PW>=Pwdiv2,1,'first');
    UnderHalfIdx(:,i)= find(PW<=Pwdiv2,1,'last');
    MedianPSD(:,i) = (Freq(OverHalfIdx(:,i))+Freq(UnderHalfIdx(:,i)))/2

[Star Strider]

I only suggested them to you.

Use whatever works to do what you want.