how to pick partial data from a waveform?

Hi

I have a wave form data from a file, where I am interested in only a part of the data. Request you to help me or suggest me how to extract only a portion of the data that I am interested.

To explain better one may be able to see orange curve it has three different stages 1. sin wave with small amplitued, 2. in between non sign wave 3. sign wave with higher apmplitude.

For reference i have attached sample file with this query.

I need to extract all the data only in the 2nd stage.

I need to extract only part of the data between the two purple vertical lines as shown below.

Answer 1

[Sriramakrishna Turaga]

@Mathieu NOE

Hi I have used signal processing tool box and started working on this code. after the specgram function I have taken

sum(sg_dBpeak) 

as a measure for my data filter and from there findingout the region of interest. This is very well working for test.txt but not test2.txt.

One more important criteria for the data filteration is I can ignore all the duration of the data where current is

Current = data(:,2);

flat. so to do this I have used some threshold value to filterout that data and then started appliying the specgram on that filteredout data.

Here is the below code.

clc
clear all
close all
% data = readmatrix('test.txt',"NumHeaderLines",1);
data = readmatrix('test2.txt',"NumHeaderLines",1);

Time = data(:,1);
Current = data(:,2);
Voltage = data(:,3);
samples = length(Time);
Fs = (samples-1)/(Time(end) - Time(1));

mn=min(Current);mx=max(Current);
% threshold=1.35/100;
threshold=2.6/100;
mnth=threshold*mn
mxth=threshold*mx
goodrows=find(Current(:)>mxth | Current(:)<mnth);
df_T=Time(goodrows,:);df_I=Current(goodrows,:);df_vs=Voltage(goodrows,:);

figure(1),
subplot(2,2,1),plot(Time,Current,'r-'),subplot(2,2,2),plot(df_T,df_I,'b');
subplot(2,2,3),plot(Time,Voltage,'r-'),subplot(2,2,4),plot(df_T,df_vs,'b');

% NB : decim = 1 will do nothing (output = input)
% decim = 5;
% if decim>1
%     Current_decim = decimate(df_I,decim);
%     Voltage_decim = decimate(df_vs,decim);
%     Fs = Fs/decim;
% end
samples = length(df_vs);
% Time_decim = (0:samples-1)*1/Fs;
% specgramdemo(Voltage,Fs)
% SPECTROGRAM STFT ANALYSIS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = 100;    % 
Overlap = 0.95;
w = hanning(NFFT);     % Hanning window / Use the HANN function to get a Hanning window which has the first and last zero-weighted samples.
% spectrogram dB scale
spectrogram_dB_scale = 80;  % dB range scale (means , the lowest displayed level is XX dB below the max level)

[sg,fsg,tsg] = specgram(df_vs,NFFT,Fs,hanning(NFFT),floor(NFFT*Overlap));  

% FFT normalisation and conversion amplitude from linear to dB (peak)
sg_dBpeak = 20*log10(abs(sg))+20*log10(2/length(fsg));     % NB : X=fft(x.*hanning(N))*4/N; % hanning only


% saturation of the dB range : 
% saturation_dB = 60;  % dB range scale (means , the lowest displayed level is XX dB below the max level)
min_disp_dB = round(max(max(sg_dBpeak))) - spectrogram_dB_scale;
sg_dBpeak(sg_dBpeak<min_disp_dB) = min_disp_dB;
kk=sum(sg_dBpeak);

% plots spectrogram
figure(2);
ind = find(fsg>1);
fsg = fsg(ind);
sg_dBpeak = sg_dBpeak(ind,:);
% imagesc(tsg,fsg,sg_dBpeak);colormap('jet');
subplot(3,1,1),plot(tsg,interp1(df_T,df_vs,tsg))
subplot(3,1,2),imagesc(tsg,fsg,sg_dBpeak);colormap('jet');
axis('xy');grid
title(['Spectrogram / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(fsg(2)-fsg(1)) ' Hz ']);
xlabel('Time (s)');ylabel('Frequency (Hz)');
subplot(3,1,3),plot(tsg,kk);

% kk=sum(sg_dBpeak);
% tsg1=find(kk>500);
tsg1=find(kk>750);
tmax=max(df_T)+0.5e-3
tmin=min(tsg(tsg1))-0.5e-3


T_final=Time(find(Time(:)>tmin & Time(:)<tmax));
I_final=Current(find(Time(:)>tmin & Time(:)<tmax));
U_final=Voltage(find(Time(:)>tmin & Time(:)<tmax));
figure(3);
subplot(2,1,1),plot(Time,Voltage,'b',T_final,U_final,'r')
subplot(2,1,2),plot(Time,Current,'b',T_final,I_final,'r')

One more thing that I have observed is with the myspecgram function (function defined by you to get out of signal proceesing tool box) the

sum(sg_dBpak)

is turning out to be all -ve values.

Request you to help me from here.

Thanks & Regards

Sriramakrishna Turaga

[Sriramakrishna Turaga]

Hi @Mathieu NOE

Finally I could do the task with the following code. But this is with all the functions related to signal processing tool box (not user defined).

clc
clear all
close all
data = readmatrix('test.txt',"NumHeaderLines",1);
% data = readmatrix('test2.txt',"NumHeaderLines",1);

Time = data(:,1);
Current = data(:,2);
Voltage = data(:,3);
samples = length(Time);
Fs = (samples-1)/(Time(end) - Time(1));
C_s=movmean(Current,100);V_s=movmean(Voltage,100);
mn=min(C_s);mx=max(C_s);
threshold=1.35/100;
% threshold=2.6/100;
mnth=threshold*mn
mxth=threshold*mx
goodrows=find(C_s(:)>mxth | C_s(:)<mnth);
% df_T=Time(goodrows,:);df_I=Current(goodrows,:);df_vs=V_s(goodrows,:);
df_T=Time(goodrows,:);df_I=Current(goodrows,:);df_vs=Voltage(goodrows,:);
figure(1),
subplot(2,2,1),plot(Time,Current,'r-'),subplot(2,2,2),plot(df_T,df_I,'b');
subplot(2,2,3),plot(Time,Voltage,'r-'),subplot(2,2,4),plot(df_T,df_vs,'b');

% NB : decim = 1 will do nothing (output = input)
decim = 1;
if decim>1
    Current_decim = decimate(df_I,decim);
    Voltage_decim = decimate(df_vs,decim);
else
    Current_decim = df_I;
    Voltage_decim = df_vs;   
end
Fs = Fs/decim;
samples = length(Voltage_decim);
Time_decim = (0:samples-1)*1/Fs;
% specgramdemo(Voltage,Fs)
% SPECTROGRAM STFT ANALYSIS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = 100;    % 
Overlap = 0.95;
w = hanning(NFFT);     % Hanning window / Use the HANN function to get a Hanning window which has the first and last zero-weighted samples.
% spectrogram dB scale
spectrogram_dB_scale = 80;  % dB range scale (means , the lowest displayed level is XX dB below the max level)

[sg,fsg,tsg] = specgram(Voltage_decim,NFFT,Fs,hanning(NFFT),floor(NFFT*Overlap));  

% FFT normalisation and conversion amplitude from linear to dB (peak)
sg_dBpeak = 20*log10(abs(sg))+20*log10(2/length(fsg));     % NB : X=fft(x.*hanning(N))*4/N; % hanning only


% saturation of the dB range : 
% saturation_dB = 60;  % dB range scale (means , the lowest displayed level is XX dB below the max level)
min_disp_dB = round(max(max(sg_dBpeak))) - spectrogram_dB_scale;
sg_dBpeak(sg_dBpeak<min_disp_dB) = min_disp_dB;
kk=sum(sg_dBpeak);


% plots spectrogram
figure(2);
ind = find(fsg>1);
fsg = fsg(ind);
sg_dBpeak = sg_dBpeak(ind,:);
% imagesc(tsg,fsg,sg_dBpeak);colormap('jet');
subplot(3,1,1),plot(Time_decim,Voltage_decim)
subplot(3,1,2),imagesc(tsg,fsg,sg_dBpeak);colormap('jet');
axis('xy');grid
title(['Spectrogram / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(fsg(2)-fsg(1)) ' Hz ']);
xlabel('Time (s)');ylabel('Frequency (Hz)');
subplot(3,1,3),plot(tsg,kk);


% kk=sum(sg_dBpeak);
% tsg1=find(kk>500);
tsg1=tsg(10:end);
kk=kk(10:end);
tsg1=find(kk>500);
tmax=max(df_T)+0.5e-3
tmin=min(tsg(kk>700))-0.5e-3
% 
% 
T_final=Time(find(Time_decim(:)>tmin & Time_decim(:)<tmax));
I_final=Current_decim(find(Time_decim(:)>tmin & Time_decim(:)<tmax));
U_final=Voltage_decim(find(Time_decim(:)>tmin & Time_decim(:)<tmax));
figure(3);
subplot(2,1,1),plot(Time_decim,Voltage_decim,'b',T_final,U_final,'r')
subplot(2,1,2),plot(Time_decim,Current_decim,'b',T_final,I_final,'r')

Please find the attached plots where required signal is extracted. Request you to help me with the same with the functions that doesn't require any toolbox or only core matlab functions specially 1. decimate, 2. Specgram.

What i have observed is the mydecimate function and myspecgram function you gave me are not yielding the same result i couldn't understand the reason all the sg_dBpeak values are -ve where as with the functions from signal processing tool box are positive, this is makeing me difficult to filter the data based on the

sum(sg_dBpeak)

function between the two (signal processing tool box vs core matlab functions).

Plots with test.txt data

Data from the file and after removing the current zero portion below

Final signal overlapped

plots with test2.txt

Data from the file and after removing the current zero portion below

Final data overlapped

[Mathieu NOE]

Glad you have now a working solution !

Good luck for the future !

Answer 2

[Mathieu NOE]

helloo

this is my submission, hope it helps !

clc
clear all
close all

data = readmatrix('test.txt',"NumHeaderLines",1);
% [data,HEAD] = readclm('test.txt',3,1);
x = data(:,1);
y = data(:,2);
z = data(:,3);

% smoothin a bit z to reduce noise in the selected area (to avoid multi^ple
% crossing points selection - below)
zs = movmean(z,50);

figure(1),
subplot(2,1,1),plot(x,y)
subplot(2,1,2),plot(x,z,'b',x,zs,'r')

threshold = 150; % user defined

[ind_pos,t0_pos,s0_pos,ind_neg,t0_neg,s0_neg]= crossing_V7(zs,x,threshold,'linear'); % positive (pos) and negative (neg) slope crossing points 
% ind => time index (samples)
% t0 => corresponding time (x) values 
% s0 => corresponding function (y) values , obviously they must be equal to "threshold"

figure(2)
plot(x,z,t0_pos,s0_pos,'+r',t0_neg,s0_neg,'+g','linewidth',2,'markersize',12);grid
legend('signal','positive slope crossing points','negative slope crossing points');

% selected portion is between first pos slope and first neg slope crossing
% points plus some pre and post buffer 
buffer = 0.0005; % in seconds
ind = find(x>=t0_pos(1)-buffer & x<=t0_neg(1)+buffer);

figure(3)
plot(x(ind),z(ind));
title('Selected / Extracted signal')


function [ind_pos,t0_pos,s0_pos,ind_neg,t0_neg,s0_neg] = crossing_V7(S,t,level,imeth)

% [ind,t0,s0,t0close,s0close] = crossing_V6(S,t,level,imeth,slope_sign) % older format

% CROSSING find the crossings of a given level of a signal
%   ind = CROSSING(S) returns an index vector ind, the signal
%   S crosses zero at ind or at between ind and ind+1
%   [ind,t0] = CROSSING(S,t) additionally returns a time
%   vector t0 of the zero crossings of the signal S. The crossing
%   times are linearly interpolated between the given times t
%   [ind,t0] = CROSSING(S,t,level) returns the crossings of the
%   given level instead of the zero crossings
%   ind = CROSSING(S,[],level) as above but without time interpolation
%   [ind,t0] = CROSSING(S,t,level,par) allows additional parameters
%   par = {'none'|'linear'}.
%	With interpolation turned off (par = 'none') this function always
%	returns the value left of the zero (the data point thats nearest
%   to the zero AND smaller than the zero crossing).
%
%	[ind,t0,s0] = ... also returns the data vector corresponding to 
%	the t0 values.
%
%	[ind,t0,s0,t0close,s0close] additionally returns the data points
%	closest to a zero crossing in the arrays t0close and s0close.
%
%	This version has been revised incorporating the good and valuable
%	bugfixes given by users on Matlabcentral. Special thanks to
%	Howard Fishman, Christian Rothleitner, Jonathan Kellogg, and
%	Zach Lewis for their input. 

% Steffen Brueckner, 2002-09-25
% Steffen Brueckner, 2007-08-27		revised version

% Copyright (c) Steffen Brueckner, 2002-2007
% [e-mail]

% M Noe
% added positive or negative slope condition

% check the number of input arguments
error(nargchk(1,4,nargin));

% check the time vector input for consistency
if nargin < 2 | isempty(t)
	% if no time vector is given, use the index vector as time
    t = 1:length(S);
elseif length(t) ~= length(S)
	% if S and t are not of the same length, throw an error
    error('t and S must be of identical length!');    
end

% check the level input
if nargin < 3
	% set standard value 0, if level is not given
    level = 0;
end

% check interpolation method input
if nargin < 4
    imeth = 'linear';
end


% make row vectors
t = t(:)';
S = S(:)';

% always search for zeros. So if we want the crossing of 
% any other threshold value "level", we subtract it from
% the values and search for zeros.
S   = S - level;

% first look for exact zeros
ind0 = find( S == 0 ); 

% then look for zero crossings between data points
S1 = S(1:end-1) .* S(2:end);
ind1 = find( S1 < 0 );

% bring exact zeros and "in-between" zeros together 
ind = sort([ind0 ind1]);

% and pick the associated time values
t0 = t(ind); 
s0 = S(ind);

if strcmp(imeth,'linear')
    % linear interpolation of crossing
    for ii=1:length(t0)
        %if abs(S(ind(ii))) > eps(S(ind(ii)))    % MATLAB V7 et +
        if abs(S(ind(ii))) > eps*abs(S(ind(ii)))    % MATLAB V6 et -    EPS * ABS(X)
            
            % interpolate only when data point is not already zero
            NUM = (t(ind(ii)+1) - t(ind(ii)));
            DEN = (S(ind(ii)+1) - S(ind(ii)));
            slope =  NUM / DEN;
            slope_sign(ii) = sign(slope);
            t0(ii) = t0(ii) - S(ind(ii)) * slope;
            s0(ii) = level;
        end
    end
end

% extract the positive slope crossing points 
ind_pos = find(sign(slope_sign)>0);
t0_pos = t0(ind_pos);
s0_pos = s0(ind_pos);

% extract the negative slope crossing points 
ind_neg = find(sign(slope_sign)<0);
t0_neg = t0(ind_neg);
s0_neg = s0(ind_neg);

end

% % Addition:
% % Some people like to get the data points closest to the zero crossing,
% % so we return these as well
% [CC,II] = min(abs([S(ind-1) ; S(ind) ; S(ind+1)]),[],1); 
% ind2 = ind + (II-2); %update indices 
% 
% t0close = t(ind2);
% s0close = S(ind2);

[Sriramakrishna Turaga]

@Mathieu NOE

Thank you very much for your support. looks like it is working for this case, but i have many such data files one more attached here. it is not able to extract similarly...

i wanted to extract all other curves in the same range (not the just one voltage curve). Request you help me understand the threshold value little more detailed.

Thanks & Regards

Sriramakrishna Turaga

[Mathieu NOE]

hello

so this is a completely different case

which section of the signal you want to extract ? is there a way to define this section of data for all your files ?

all the best

Answer 3

[Sriramakrishna Turaga]

Hi @Mathieu NOE

Yes this a different case. The way we define section of data is to identify a location it is not in sine wave form.

All my wave forms will be one way similar, where the orange curve (voltage plot) will have non-sinusudol form as depicted in the picture attached. in the previous picutre also I want to extract only that data where the voltage deviates from its sine wave form.

So the challenge is to identify from which point of time or instance the curve is deviating from being sine wave form.

Thanks in advance.

Sriramakrishna Turaga

[Mathieu NOE]

ok , I understand better

there is a difficulty when we think only in terms of sinus / not sinus , is that there are also portion of signals where transients are decaying waves with higher frequency sinus content

I'm moving now using STFT / spectrogram to defne the sinus / non sinus areas and it maybe the way to go

[Mathieu NOE]

now this needs to me coded :

Answer 4

[Sriramakrishna Turaga]

oh great!!!

When I am going through the matlab functions available (I am very much newbi to matlab). For this kind of data extraction I am thinking of dividing the data in to small small chunks (one chunk is one cycle of data). perform fft supress all the sin waves take inverse fft after suppressing the sin waves. That would be almost close to what I am looking for...

But this method seems to be more visual and elegent unfortunately I don't have signal processing tools box so what ever needs to be done i need to stick to core matlab functions only. I think you gave me the idea that is almost close to what I am looking for.

I would like to explain little more insite on the type of data that I get.

1. the first one test.txt file we call it as O-type. In this the voltage wave initially will be of small amplitude and directly goes into the data region what I am interested in.
2. this is one more type of data test2.txt file and CO-type. In this the voltage initially with big amplitude, after few cycles experiment starts where for about 100 milli seconds voltage will be a second wave form, after that only it enters into the data region that we are looking for.

And in all the cases that data to be extracted will be spread across min 10-15 milli seconds of time span.

I think with an extra input from the user on the type of data O-type or CO-type and non-sinus region that is spread across more than 10-15 ms will surely give us the data to be extracted.

Request you to guide me for the same.

Thank you very much for your support

Sriramakrishna Turaga

[Mathieu NOE]

hello - I'll try to help yopu the best I can

yes we can probably get the best results if we can tell the code what kind of signal is feeded and what should be the min / max duration of the extracted signal (should it be only one segment ?)

FYI, this is the code below , but I can easily rewritte the spectrogram function using the fft function :

data = readmatrix('test.txt',"NumHeaderLines",1);

% data = readmatrix('test2.txt',"NumHeaderLines",1);
% [data,HEAD] = readclm('test.txt',3,1);
time = data(:,1);
y = data(:,2);
z = data(:,3);
samples = length(time);
Fs = (samples-1)/(time(end) - time(1));

figure(1),
subplot(2,1,1),plot(time,y)
subplot(2,1,2),plot(time,z)

% NB : decim = 1 will do nothing (output = input)
decim = 10;
if decim>1
    y = decimate(y,decim);
    z = decimate(z,decim);
    Fs = Fs/decim;
end
samples = length(z);
time = (0:samples-1)*1/Fs;
% specgramdemo(z,Fs)
% SPECTROGRAM STFT ANALYSIS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = 100;    % 
Overlap = 0.95;
w = hanning(NFFT);     % Hanning window / Use the HANN function to get a Hanning window which has the first and last zero-weighted samples.
% spectrogram dB scale
spectrogram_dB_scale = 80;  % dB range scale (means , the lowest displayed level is XX dB below the max level)

[sg,fsg,tsg] = specgram(z,NFFT,Fs,hanning(NFFT),floor(NFFT*Overlap));  

% FFT normalisation and conversion amplitude from linear to dB (peak)
sg_dBpeak = 20*log10(abs(sg))+20*log10(2/length(fsg));     % NB : X=fft(x.*hanning(N))*4/N; % hanning only


% saturation of the dB range : 
% saturation_dB = 60;  % dB range scale (means , the lowest displayed level is XX dB below the max level)
min_disp_dB = round(max(max(sg_dBpeak))) - spectrogram_dB_scale;
sg_dBpeak(sg_dBpeak<min_disp_dB) = min_disp_dB;

% plots spectrogram
figure(2);
ind = find(fsg>1);
fsg = fsg(ind);
sg_dBpeak = sg_dBpeak(ind,:);
% imagesc(tsg,fsg,sg_dBpeak);colormap('jet');
subplot(2,1,1),plot(tsg,interp1(time,z,tsg))
subplot(2,1,2),imagesc(tsg,fsg,sg_dBpeak);colormap('jet');
axis('xy');grid
title(['Spectrogram / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(fsg(2)-fsg(1)) ' Hz ']);
xlabel('Time (s)');ylabel('Frequency (Hz)');

[Mathieu NOE]

this is a bit improved code , and also I wrote my own spectrogram function based on fft

clc
clear all
close all

data = readmatrix('test.txt',"NumHeaderLines",1);
% data = readmatrix('test2.txt',"NumHeaderLines",1);

time = data(:,1);
y = data(:,2);
z = data(:,3);
samples = length(time);
Fs = (samples-1)/(time(end) - time(1));

% figure(1),
% subplot(2,1,1),plot(time,y)
% subplot(2,1,2),plot(time,z)
% 
% NB : decim = 1 will do nothing (output = input)
decim = 5;
if decim>1
    y = decimate(y,decim);
    z = decimate(z,decim);
    Fs = Fs/decim;
end
samples = length(z);
time = (0:samples-1)*1/Fs;

% SPECTROGRAM STFT ANALYSIS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = 128;    % 
Overlap = 0.95;
w = hanning(NFFT);     % Hanning window / Use the HANN function to get a Hanning window which has the first and last zero-weighted samples.
% spectrogram dB scale
spectrogram_dB_scale = 80;  % dB range scale (means , the lowest displayed level is XX dB below the max level)

% [sg,fsg,tsg] = specgram(z,NFFT,Fs,hanning(NFFT),floor(NFFT*Overlap));  
[sg,fsg,tsg]  = myspecgram(z,Fs,NFFT,Overlap);  % home made spectrogram function

% FFT normalisation and conversion amplitude from linear to dB (peak)
sg_dBpeak = 20*log10(abs(sg))+20*log10(2/length(fsg));     % NB : X=fft(x.*hanning(N))*4/N; % hanning only


% saturation of the dB range : 
% saturation_dB = 60;  % dB range scale (means , the lowest displayed level is XX dB below the max level)
min_disp_dB = round(max(max(sg_dBpeak))) - spectrogram_dB_scale;
sg_dBpeak(sg_dBpeak<min_disp_dB) = min_disp_dB;

% plots spectrogram
figure(2);
ind = find(fsg>1);
fsg = fsg(ind);
sg_dBpeak = sg_dBpeak(ind,:);
subplot(2,1,1),plot(time,z)
xlim([min(tsg) max(tsg)]);
subplot(2,1,2),imagesc(tsg,fsg,sg_dBpeak);colormap('jet');
xlim([min(tsg) max(tsg)]);
axis('xy');grid
title(['Spectrogram / Fs = ' num2str(round(Fs)) ' Hz / Delta f = ' num2str(round(fsg(2)-fsg(1))) ' Hz ']);
xlabel('Time (s)');ylabel('Frequency (Hz)');


function  [fft_specgram,freq_vector,time] = myspecgram(signal, Fs, nfft, Overlap)
% FFT peak spectrogram of signal  (example sinus amplitude 1   = 0 dB after fft).
%   signal - input signal, 
%   Fs - Sampling frequency (Hz).
%   nfft - FFT window size
%   Overlap - buffer overlap % (between 0 and 0.95)
signal = signal(:);
samples = length(signal);
% fill signal with zeros if its length is lower than nfft
if samples<nfft
    s_tmp = zeros(nfft,1);
    s_tmp((1:samples)) = signal;
    signal = s_tmp;
    samples = nfft;
end

% window : hanning
window = hanning(nfft);
window = window(:);

%    compute fft with overlap 
 offset = fix((1-Overlap)*nfft);
 spectnum = 1+ fix((samples-nfft)/offset); % Number of windows
%     % for info is equivalent to : 
%     noverlap = Overlap*nfft;
%     spectnum = fix((samples-noverlap)/(nfft-noverlap));	% Number of windows

    % main loop
    fft_specgram = [];
    for ci=1:spectnum
        start = (ci-1)*offset;
        sw = signal((1+start):(start+nfft)).*window;
        fft_specgram = [fft_specgram abs(fft(sw))*4/nfft];     % X=fft(x.*hanning(N))*4/N; % hanning only 
    end
    fft_specgram = fft_specgram/spectnum; % to do linear averaging scaling

% one sidded fft spectrum  % Select first half 
    if rem(nfft,2)    % nfft odd
        select = (1:(nfft+1)/2)';
    else
        select = (1:nfft/2+1)';
    end
fft_specgram = fft_specgram(select,:);
freq_vector = (select - 1)*Fs/nfft;

% time vector 
% time stamps are defined in the middle of the buffer
time = ((0:spectnum-1)*offset + round(offset/2))/Fs;

end

[Sriramakrishna Turaga]

Hi @Mathieu NOE

Thank you very much for your support, in both the versions the decimate is a signal processing tool box specific function, hence I am not able to test. So when you tried at your end this code is able to extract the required data only?

[Mathieu NOE]

this is the code with correction for missing decimate function

but no I am not yet to the point of extracting the data

clc
clear all
close all

data = readmatrix('test.txt',"NumHeaderLines",1);
% data = readmatrix('test2.txt',"NumHeaderLines",1);

time = data(:,1);
y = data(:,2);
z = data(:,3);
samples = length(time);
Fs = (samples-1)/(time(end) - time(1));

% figure(1),
% subplot(2,1,1),plot(time,y)
% subplot(2,1,2),plot(time,z)
% 
% NB : decim = 1 will do nothing (output = input)
decim = 5;
if decim>1
    y = mydecimate(y,decim);
    z = mydecimate(z,decim);
    Fs = Fs/decim;
end
samples = length(z);
time = (0:samples-1)*1/Fs;

% SPECTROGRAM STFT ANALYSIS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = 128;    % 
Overlap = 0.95;
w = hanning(NFFT);     % Hanning window / Use the HANN function to get a Hanning window which has the first and last zero-weighted samples.
% spectrogram dB scale
spectrogram_dB_scale = 80;  % dB range scale (means , the lowest displayed level is XX dB below the max level)

% [sg,fsg,tsg] = specgram(z,NFFT,Fs,hanning(NFFT),floor(NFFT*Overlap));  
[sg,fsg,tsg]  = myspecgram(z,Fs,NFFT,Overlap);  % home made spectrogram function

% FFT normalisation and conversion amplitude from linear to dB (peak)
sg_dBpeak = 20*log10(abs(sg))+20*log10(2/length(fsg));     % NB : X=fft(x.*hanning(N))*4/N; % hanning only


% saturation of the dB range : 
% saturation_dB = 60;  % dB range scale (means , the lowest displayed level is XX dB below the max level)
min_disp_dB = round(max(max(sg_dBpeak))) - spectrogram_dB_scale;
sg_dBpeak(sg_dBpeak<min_disp_dB) = min_disp_dB;

% plots spectrogram
figure(2);
ind = find(fsg>1);
fsg = fsg(ind);
sg_dBpeak = sg_dBpeak(ind,:);
subplot(2,1,1),plot(time,z)
xlim([min(tsg) max(tsg)]);
subplot(2,1,2),imagesc(tsg,fsg,sg_dBpeak);colormap('jet');
xlim([min(tsg) max(tsg)]);
axis('xy');grid
title(['Spectrogram / Fs = ' num2str(round(Fs)) ' Hz / Delta f = ' num2str(round(fsg(2)-fsg(1))) ' Hz ']);
xlabel('Time (s)');ylabel('Frequency (Hz)');


function  [fft_specgram,freq_vector,time] = myspecgram(signal, Fs, nfft, Overlap)
% FFT peak spectrogram of signal  (example sinus amplitude 1   = 0 dB after fft).
%   signal - input signal, 
%   Fs - Sampling frequency (Hz).
%   nfft - FFT window size
%   Overlap - buffer overlap % (between 0 and 0.95)
signal = signal(:);
samples = length(signal);
% fill signal with zeros if its length is lower than nfft
if samples<nfft
    s_tmp = zeros(nfft,1);
    s_tmp((1:samples)) = signal;
    signal = s_tmp;
    samples = nfft;
end

% window : hanning
window = hanning(nfft);
window = window(:);

%    compute fft with overlap 
 offset = fix((1-Overlap)*nfft);
 spectnum = 1+ fix((samples-nfft)/offset); % Number of windows
%     % for info is equivalent to : 
%     noverlap = Overlap*nfft;
%     spectnum = fix((samples-noverlap)/(nfft-noverlap));	% Number of windows

    % main loop
    fft_specgram = [];
    for ci=1:spectnum
        start = (ci-1)*offset;
        sw = signal((1+start):(start+nfft)).*window;
        fft_specgram = [fft_specgram abs(fft(sw))*4/nfft];     % X=fft(x.*hanning(N))*4/N; % hanning only 
    end
    fft_specgram = fft_specgram/spectnum; % to do linear averaging scaling

% one sidded fft spectrum  % Select first half 
    if rem(nfft,2)    % nfft odd
        select = (1:(nfft+1)/2)';
    else
        select = (1:nfft/2+1)';
    end
fft_specgram = fft_specgram(select,:);
freq_vector = (select - 1)*Fs/nfft;

% time vector 
% time stamps are defined in the middle of the buffer
time = ((0:spectnum-1)*offset + round(offset/2))/Fs;

end


function odata = mydecimate(idata,r,nfilt,option)
%DECIMATE Resample data at a lower rate after lowpass filtering.
%   Y = DECIMATE(X,R) resamples the sequence in vector X at 1/R times the
%   original sample rate.  The resulting resampled vector Y is R times
%   shorter, i.e., LENGTH(Y) = CEIL(LENGTH(X)/R). By default, DECIMATE
%   filters the data with an 8th order Chebyshev Type I lowpass filter with
%   cutoff frequency .8*(Fs/2)/R, before resampling.
%
%   Y = DECIMATE(X,R,N) uses an N'th order Chebyshev filter.  For N greater
%   than 13, DECIMATE will produce a warning regarding the unreliability of
%   the results.  See NOTE below.
%
%   Y = DECIMATE(X,R,'FIR') uses a 30th order FIR filter generated by
%   FIR1(30,1/R) to filter the data.
%
%   Y = DECIMATE(X,R,N,'FIR') uses an Nth FIR filter.
%
%   Note: For better results when R is large (i.e., R > 13), it is
%   recommended to break R up into its factors and calling DECIMATE several
%   times.
%
%   EXAMPLE: Decimate a signal by a factor of four
%   t = 0:.00025:1;  % Time vector
%   x = sin(2*pi*30*t) + sin(2*pi*60*t);
%   y = decimate(x,4);
%   subplot(1,2,1);
%   stem(x(1:120)), axis([0 120 -2 2])   % Original signal
%   title('Original Signal')
%   subplot(1,2,2);
%   stem(y(1:30))                        % Decimated signal
%   title('Decimated Signal')
%
%   See also DOWNSAMPLE, INTERP, RESAMPLE, FILTFILT, FIR1, CHEBY1.

%   Author(s): L. Shure, 6-4-87
%              L. Shure, 12-11-88, revised
%   Copyright 1988-2018 The MathWorks, Inc.

%   References:
%    [1] "Programs for Digital Signal Processing", IEEE Press
%         John Wiley & Sons, 1979, Chap. 8.3.

narginchk(2,4);
error(nargoutchk(0,1,nargout,'struct'));

if nargin > 2
    nfilt = convertStringsToChars(nfilt);
end

if nargin > 3
    option = convertStringsToChars(option);
end

% Validate required inputs 
validateinput(idata,r);

if fix(r) == 1
    odata = idata;
    return
end

fopt = 1;
if nargin == 2
    nfilt = 8;
else
    if ischar(nfilt)
        if nfilt(1) == 'f' || nfilt(1) == 'F'
            fopt = 0;
        elseif nfilt(1) == 'i' || nfilt(1) == 'I'
            fopt = 1;
        else
            error(message('signal:decimate:InvalidEnum'))
        end
        if nargin == 4
            nfilt = option;
        else
            nfilt = 8*fopt + 30*(1-fopt);
        end
    else
        if nargin == 4
            if option(1) == 'f' || option(1) == 'F'
                fopt = 0;
            elseif option(1) == 'i' || option(1) == 'I'
                fopt = 1;
            else
                error(message('signal:decimate:InvalidEnum'))
            end
        end
    end
end
if fopt == 1 && nfilt > 13
    warning(message('signal:decimate:highorderIIRs'));
end

nd = length(idata);
m = size(idata,1);
nout = ceil(nd/r);

if fopt == 0	% FIR filter
    b = fir1(nfilt,1/r);
    % prepend sequence with mirror image of first points so that transients
    % can be eliminated. then do same thing at right end of data sequence.
    nfilt = nfilt+1;
    itemp = 2*idata(1) - idata((nfilt+1):-1:2);
    [itemp,zi]=filter(b,1,itemp); %#ok
    [odata,zf] = filter(b,1,idata,zi);
    if m == 1	% row data
        itemp = zeros(1,2*nfilt);
    else	% column data
        itemp = zeros(2*nfilt,1);
    end
    itemp(:) = 2*idata(nd)-idata((nd-1):-1:(nd-2*nfilt));
    itemp = filter(b,1,itemp,zf);
    % finally, select only every r'th point from the interior of the lowpass
    % filtered sequence
    gd = grpdelay(b,1,8);
    list = round(gd(1)+1.25):r:nd;
    odata = odata(list);
    lod = length(odata);
    nlen = nout - lod;
    nbeg = r - (nd - list(length(list)));
    if m == 1	% row data
        odata = [odata itemp(nbeg:r:nbeg+nlen*r-1)];
    else
        odata = [odata; itemp(nbeg:r:nbeg+nlen*r-1)];
    end
else	% IIR filter
    rip = .05;	% passband ripple in dB
    [b,a] = cheby1(nfilt, rip, .8/r);
    while all(b==0) || (abs(filtmag_db(b,a,.8/r)+rip)>1e-6)
        nfilt = nfilt - 1;
        if nfilt == 0
            break
        end
        [b,a] = cheby1(nfilt, rip, .8/r);
    end
    if nfilt == 0
        error(message('signal:decimate:InvalidRange'))
    end

    % be sure to filter in both directions to make sure the filtered data has zero phase
    % make a data vector properly pre- and ap- pended to filter forwards and back
    % so end effects can be obliterated.
    odata = filtfilt(b,a,idata);
    nbeg = r - (r*nout - nd);
    odata = odata(nbeg:r:nd);
end

%--------------------------------------------------------------------------
function H = filtmag_db(b,a,f)
%FILTMAG_DB Find filter's magnitude response in decibels at given frequency.

nb = length(b);
na = length(a);
top = exp(-1i*(0:nb-1)*pi*f)*b(:);
bot = exp(-1i*(0:na-1)*pi*f)*a(:);

H = 20*log10(abs(top/bot));
end
%--------------------------------------------------------------------------
function validateinput(x,r)
% Validate 1st two input args: signal and decimation factor

if isempty(x) || issparse(x) || ~isa(x,'double')
    error(message('signal:decimate:invalidInput', 'X'));
end

if (abs(r-fix(r)) > eps) || (r <= 0)
    error(message('signal:decimate:invalidR', 'R'));
end
end
end