How to load a specific column from files + FFT function

Hi, I would like to write a script that reads only the 7th column from text files (columns are separated by semicolon ";" ). Additionally, I need to count the FFT of each 7th column. I have started to do this but something is not working for me. Maybe someone has a better idea for this code.

%%%%%
dane = dir('srednie'); %srednie is the name of the folder 
dane = dane(3:8);

for k = 1:1:max(size(dane))
    
    file = horzcat('srednie/', dane(k).name);
    [id, kom] = fopen(file);
    if id <0
        disp(kom)
        
    end
    
    MyData.(horzcat('srednie_', dane(k).name(1:end-4))) = ...
        fscanf(id,'%f;%f;%f;%f;%f;%f;%f;%f', [8 7200]);
end


dataStruct = ????;  %how to load files from a structure?
                                   
result = zeros(6,7200);
dataNames = fieldnames(dataStruct);
    for i = 1:length(dataNames)
        result{i,1} = dataStruct.dataNames{i}(7,:);
    end

[Michael]

Thank you for your answer!

Answer 1

[Mathieu NOE]

helo again

as the signals are quite nn stationnary I assumed that a spectrogram is more appropriate than a "simple" FFT (with or without averaging)

try this and let me know what you think

the spectrogram amplitude are converted into dB (log rangeà , otherwise in linear range it's rather difficult to distinguish the small amplitudes beside the large ones.

the x axis is the angle 0 to 720° (could be also displayed as time increment, this corresponds to the commented code FYI)

all the best

% These are data concerning engine vibrations.
% One file is one cycle. In one cycle the shaft rotates twice,
% i.e. 720 degrees, the measurement was made 7200 times per cycle,
% i.e. every 0.1 degree of rotation. 

% FFT / spectrogram  parameters
RPM = 2000; % my assumption
Fs = RPM/60*3600; % 3600 samples for one shaft rev
dt = 1/Fs;
nfft = 500; % frame length
overlap = 0.75; % window overlap % overlap is 95% of nfft here

% read current filenames  in folder
S = dir('Z_data*.txt');
S = natsortfiles(S); % natsortfiles now works on the entire structure
% natsortfiles available from FEX : 
% https://fr.mathworks.com/matlabcentral/fileexchange/47434-natural-order-filename-sort?s_tid=srchtitle
figure(1),hold on 
for k = 1:numel(S)
    S(k).name % display file name in command window (for info only, you can remove / comment this line)
    F = fullfile(S(k).folder, S(k).name);
    out = readmatrix(F);
    [samples,cols] = size(out);
    angl = 0.1*(0:samples-1); % see above : One file is one cycle. In one cycle the shaft rotates twice,
                        % i.e. 720 degrees, the measurement was made 7200 times per cycle,
                        % i.e. every 0.1 degree of rotation.
    time = dt*(0:samples-1);
    %    There are two types of files. One type has 8 columns (I am interested in the 7th column) 
%   and the other type has 4 columns (I am interested in the 4th column). 

    if cols>4 % 8 columns file 
        out = out(:,7); % this store the 7th column
    else % 4 columns file
        out = out(:,4); % this store the 4th column
    end
    
    % FFT / spectrogram analysis
    [Sout,F,T] = myspecgram(out, Fs, nfft, overlap); 
    figure(k), 
%     subplot(2,1,1),plot(time,out); % if you want the x axis labelled in time units
%     xlabel('Time (secs)') 
    subplot(2,1,1),plot(angl,out); % if you want the x axis labelled in angular units
    xlabel('Angle (°)') 
    ylabel('Amplitude')     
    xlim([0 max(angl)]);
%     subplot(2,1,2);imagesc(T,F/1000,20*log10(abs(Sout))) % dB scale / x axis labelled in time units
%     set(gca,'YDir','Normal') 
%     xlim([0 max(time)]);
%     colorbar('East');
%     xlabel('Time (secs)') 
    angl_spec = 0.1*T/dt; % create specific angle x data from time stamps of spectrogram function
    subplot(2,1,2);imagesc(angl_spec,F/1000,20*log10(abs(Sout))) % dB scale / x axis labelled in time units
    set(gca,'YDir','Normal') 
    xlim([0 max(angl)]);
    colorbar('East');
    xlabel('Angle (°)') 
    ylabel('Freq (kHz)') 
    title('Short-time Fourier Transform spectrum (dB Scale)')
    colormap('jet');
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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)


% make signal col vector
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);%    compute fft with overlap 
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

% 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(nfft/2))/Fs;

end

Answer 2

[Mathieu NOE]

hello

my first suggestion is for looping inside a folder to load txt file data

% read current filenames  in folder
S = dir('**/*.txt');

S = natsortfiles(S); % natsortfiles now works on the entire structure
% natsortfiles available from FEX : 
% https://fr.mathworks.com/matlabcentral/fileexchange/47434-natural-order-filename-sort?s_tid=srchtitle

figure(1),hold on 

for k = 1:numel(S)
    S(k).name % display file name in command window (for info only, you can remove / comment this line)
    F = fullfile(S(k).folder, S(k).name);
    %S(k).data = load(F); % or READTABLE or whatever.
    out = readmatrix(F ,"NumHeaderLines", 1);
    S(k).data = out(:,13); % this store the 13th column
    
    % plot (for fun)
    legstr{k} = S(k).name; % legend string
    plot(S(k).data);
end
legend(legstr);

% % Take a look in the structure S: it contains all of your file data and the corresponding filenames, just as you require.
% % For example, the 2nd filename and its data:
% S(2).name
% S(2).data

my second suggestion is for doing the spectral analysis of data

clc

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% data
% data = readmatrix('data2.txt',"NumHeaderLines",1);
% data = readmatrix('inputsig_non-uniform.txt');
data = readmatrix('TEC.txt');

time = data(:,1);
signal = data(:,2);

time = time(:);  % column oriented 
signal = signal(:);  % column oriented

samples = length(signal);
Fs = (samples-1)/(time(samples)-time(1)); % sampling frequency (Hz)

%% decimate (if needed)
% NB : decim = 1 will do nothing (output = input)
decim = 1;
if decim>1
    signal = decimate(signal,decim);
    Fs = Fs/decim;
end
samples = length(signal);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = 1024;    % 
OVERLAP = 0.95;

% spectrogram dB scale
spectrogram_dB_scale = 80;  % dB range scale (means , the lowest displayed level is XX dB below the max level)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% options 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% if you are dealing with acoustics, you may wish to have A weighted
% spectrums 

% option_w = 0 : linear spectrum (no weighting dB (L) )
% option_w = 1 : A weighted spectrum (dB (A) )
option_w = 0;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 1 : averaged FFT spectrum
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[freq, sensor_spectrum] = myfft_peak(signal,Fs,NFFT,OVERLAP);

% convert to dB scale (ref = 1)
sensor_spectrum_dB = 20*log10(sensor_spectrum);

% apply A weigthing if needed
if option_w == 1
    pondA_dB = pondA_function(freq);
    sensor_spectrum_dB = sensor_spectrum_dB+pondA_dB;
    my_ylabel = ('Amplitude (dB (A))');
else
    my_ylabel = ('Amplitude (dB (L))');
end


figure(1),semilogx(freq,sensor_spectrum_dB,'b');grid
title(['Averaged FFT Spectrum  / Fs = ' num2str(0.1*round(10*Fs)) ' Hz / Delta f = ' num2str(0.1*round(10*(freq(2)-freq(1)))) ' Hz ']);
xlabel('Frequency (Hz)');ylabel(my_ylabel);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 2 : time / frequency analysis : spectrogram demo
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[sg,fsg,tsg] = specgram(signal,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

% apply A weigthing if needed
if option_w == 1
    pondA_dB = pondA_function(fsg);
    sg_dBpeak = sg_dBpeak+(pondA_dB*ones(1,size(sg_dBpeak,2)));
    my_title = ('Spectrogram (dB (A))');
else
    my_title = ('Spectrogram (dB (L))');
end

% 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);
imagesc(tsg,-log10(fsg+1e-4),sg_dBpeak);colormap('jet');
axis('xy');colorbar('vert');grid
title([my_title ' / Fs = ' num2str(0.1*round(10*Fs)) ' Hz / Delta f = ' num2str(0.1*round(10*(fsg(2)-fsg(1)))) ' Hz ']);
xlabel('Time (s)');ylabel('LOG10 Period (days) ');




function pondA_dB = pondA_function(f)
	% dB (A) weighting curve
	n = ((12200^2*f.^4)./((f.^2+20.6^2).*(f.^2+12200^2).*sqrt(f.^2+107.7^2).*sqrt(f.^2+737.9^2)));
	r = ((12200^2*1000.^4)./((1000.^2+20.6^2).*(1000.^2+12200^2).*sqrt(1000.^2+107.7^2).*sqrt(1000.^2+737.9^2))) * ones(size(f));
	pondA = n./r;
	pondA_dB = 20*log10(pondA(:));
end


function  [freq_vector,fft_spectrum] = myfft_peak(signal, Fs, nfft, Overlap)
% FFT peak spectrum of signal  (example sinus amplitude 1   = 0 dB after fft).
% Linear averaging
%   signal - input signal, 
%   Fs - Sampling frequency (Hz).
%   nfft - FFT window size
%   Overlap - buffer percentage of overlap % (between 0 and 0.95)

[samples,channels] = size(signal);

% fill signal with zeros if its length is lower than nfft
if samples<nfft
    s_tmp = zeros(nfft,channels);
    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_spectrum = 0;
    for i=1:spectnum
        start = (i-1)*offset;
        sw = signal((1+start):(start+nfft),:).*(window*ones(1,channels));
        fft_spectrum = fft_spectrum + (abs(fft(sw))*4/nfft);     % X=fft(x.*hanning(N))*4/N; % hanning only 
    end
    fft_spectrum = fft_spectrum/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_spectrum = fft_spectrum(select,:);
freq_vector = (select - 1)*Fs/nfft;
end

there is not much work left to do to put the fft analysis in the first code !

[Michael]

Thank you very much for your detailed answer, I will try to apply your suggestions.

[Michael]

I have a problem, when I inserted an extra "for" loop then in the "out" there are only 4 rows and not with the necessary data. However, in the structure (variable S) there are all files and the lines are correct. Where am I making a mistake? Sorry but I am a beginner in matlab.

S = dir('srednie/*.txt');
S = natsortfiles(S); % natsortfiles now works on the entire structure 
% https://fr.mathworks.com/matlabcentral/fileexchange/47434-natural-order-filename-sort?s_tid=srchtitle
figure(1),hold off 
for k = 1:6   
    %S(k).name                                 
    F = fullfile(S(k).folder, S(k).name);
    S(k).data = load(F);                       
    out = readmatrix(F ,"NumHeaderLines", 0);  
    S(k).data = out(:,7); 
    
     for  k = 7:numel(S)
          %S(k).name                                
          F = fullfile(S(k).folder, S(k).name);
          S(k).data = load(F);                       
          out = readmatrix(F ,"NumHeaderLines", 0);  
    	  S(k).data = out(:,4);
      end    
    
   
end

[Mathieu NOE]

hello Michael

no problem, we were all beginners one day....

I suspect you don't need to have a second for loop - as far as I remember you wanted only to get the 7th column of your data ?

maybe if you zip some txt files and share them, I can better help you

all the best

[Michael]

Thank you:)

There are two types of files. One type has 8 columns (I am interested in the 7th column) and the other type has 4 columns (I am interested in the 4th column).

In the attachment there are two sample files from which I need to import data.

Z_data_8.txt - has 8th columns

Z_data_4.thx - has 4th columns

[Mathieu NOE]

hello

I am missing your sample rate to implememnt the fft analysis

FYI, this is the code to extract the data from your two files types

% read current filenames  in folder
S = dir('Z_data*.txt');
S = natsortfiles(S); % natsortfiles now works on the entire structure
% natsortfiles available from FEX : 
% https://fr.mathworks.com/matlabcentral/fileexchange/47434-natural-order-filename-sort?s_tid=srchtitle
figure(1),hold on 
for k = 1:numel(S)
    S(k).name % display file name in command window (for info only, you can remove / comment this line)
    F = fullfile(S(k).folder, S(k).name);
    out = readmatrix(F);
    [samples,cols] = size(out);

    %    There are two types of files. One type has 8 columns (I am interested in the 7th column) 
%   and the other type has 4 columns (I am interested in the 4th column). 

    if cols>4 % 8 columns file 
        out = out(:,7); % this store the 7th column
    else % 4 columns file
        out = out(:,4); % this store the 4th column
    end
    figure, plot(out)
end

[Michael]

Hello,

Thenks:) now it works!

The probing frequency was 0.1 crankshaft revolutions

[Mathieu NOE]

hello

sorry I don't understand

is that a fixed or variable acquisition rate ? is the number of samples fixed per shaft revolution ?

[Michael]

These are data concerning engine vibrations. One file is one cycle. In one cycle the shaft rotates twice, i.e. 720 degrees, the measurement was made 7200 times per cycle, i.e. every 0.1 degree of rotation.

[Mathieu NOE]

hello

I had to make some hypothesis about the sampling rate , which depends at which RPM rotates your engine

see my answer below in the "answer" section

[Michael]

WOW looks great! I must say you are good at matlab:)

I would additionally like to show the range from 300 - 500 degrees, because there is the combustion (and that is what interests me mainly).

I am also wondering how would it look like if the graph shows only maximum amplitudes in this range (as points).

Thank you very much for taking the time to tele me, I am very grateful!

[Mathieu NOE]

hello Michael

this is a first update , you can select the angle range and also the frequency range for the spectrogram

angle_range = [300 500]; % range from 300 - 500 degrees, because there is the combustion 
freq_range = [0 10]; % (units : kHz !) range to display on Frequency axis (Spectrogram y axis)

as you can see most of the energy is below 2 kHz

"I am also wondering how would it look like if the graph shows only maximum amplitudes in this range (as points)." what do you mean ? you want something like findpeaks ? (displaying the local maxima / minima ?)

full code :

% These are data concerning engine vibrations.
% One file is one cycle. In one cycle the shaft rotates twice,
% i.e. 720 degrees, the measurement was made 7200 times per cycle,
% i.e. every 0.1 degree of rotation. 

% FFT / spectrogram  parameters
RPM = 2000; % my assumption
Fs = RPM/60*3600; % 3600 samples for one shaft rev
dt = 1/Fs;
nfft = 500; % frame length
overlap = 0.75; % window overlap % overlap is 95% of nfft here
angle_range = [300 500]; % range from 300 - 500 degrees, because there is the combustion 
freq_range = [0 10]; % (units : kHz !) range to display on Frequency axis (Spectrogram y axis)

% read current filenames  in folder
S = dir('Z_data*.txt');
S = natsortfiles(S); % natsortfiles now works on the entire structure
% natsortfiles available from FEX : 
% https://fr.mathworks.com/matlabcentral/fileexchange/47434-natural-order-filename-sort?s_tid=srchtitle
figure(1),hold on 
for k = 1:numel(S)
    S(k).name % display file name in command window (for info only, you can remove / comment this line)
    F = fullfile(S(k).folder, S(k).name);
    out = readmatrix(F);
    [samples,cols] = size(out);
    angl = 0.1*(0:samples-1); % see above : One file is one cycle. In one cycle the shaft rotates twice,
    % i.e. 720 degrees, the measurement was made 7200 times per cycle,
    % i.e. every 0.1 degree of rotation.
    time = dt*(0:samples-1);
    
    %    There are two types of files. One type has 8 columns (I am interested in the 7th column) 
    %   and the other type has 4 columns (I am interested in the 4th column). 
    
    if cols>4 % 8 columns file 
        out = out(:,7); % this store the 7th column
    else % 4 columns file
        out = out(:,4); % this store the 4th column
    end
    
    % restrict to range : range from 300 - 500 degrees, because there is the combustion 
    ind = (angl>= min(angle_range) & angl<= max(angle_range)); % logical array
    angl = angl(ind);
    time = time(ind);
    out = out(ind);
    
    % FFT / spectrogram analysis
    [Sout,F,T] = myspecgram(out, Fs, nfft, overlap); 
    figure(k), 
    %     subplot(2,1,1),plot(time,out); % if you want the x axis labelled in time units
    %     xlabel('Time (secs)') 
    subplot(2,1,1),plot(angl,out); % if you want the x axis labelled in angular units
    xlabel('Angle (°)') 
    ylabel('Amplitude')     
    xlim([min(angl) max(angl)]);
    %     subplot(2,1,2);imagesc(T,F/1000,20*log10(abs(Sout))) % dB scale / x axis labelled in time units
    %     set(gca,'YDir','Normal') 
    %     xlim([0 max(time)]);
    %     colorbar('East');
    %     xlabel('Time (secs)') 
    angl_spec = (0.1*T)/dt+angl(1); % create specific angle x data from time stamps of spectrogram function
    subplot(2,1,2);imagesc(angl_spec,F/1000,20*log10(abs(Sout))) % dB scale / x axis labelled in time units
    set(gca,'YDir','Normal') 
    xlim([min(angl) max(angl)]);
    ylim([min(freq_range) max(freq_range)]);
    colorbar('East');
    xlabel('Angle (°)') 
    ylabel('Freq (kHz)') 
    title('Short-time Fourier Transform spectrum (dB Scale)')
    colormap('jet');
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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)


% make signal col vector
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);%    compute fft with overlap 
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

% 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(nfft/2))/Fs;

end

[Michael]

Thank you once again for your help. I will try to do the rest of the project myself:)

[Mathieu NOE]

OK

all the best

Answer 3

[Michael]

Hi

I have another problem. I have a series of files. There are 100 cycles stored in each file (one by one, 1 cycle to 7200, 2 cycle to 14400 etc). I need to create a matrix where:

- in the first column is the cycle number,

- in the second column max amplitude (for each cycle),

- in the third column the rotation angle at which the max amplitude occurs.

How should I do it?

[Mathieu NOE]

hello

have you started some code yet ? where have you difficulties ?

all the best

[Michael]

I do not know how I should do this. Is it better to save all cycles in separate files at the beginning and only register max amplitudes?

[Mathieu NOE]

I will try something - if you send some data files (zipped) that would be great

but I will be absent a few days during the week so I make no promise when I can answer....

[Michael]

I understand and thank you for all your help. I will try to do something myself tomorrow.

[Mathieu NOE]

hello

I thought you would send me multiples files in one zip ?

in the attached file I see 4 and not 3 columns and I am not sure to understand what I read vs. what you described above :

- in the first column is the cycle number,

- in the second column max amplitude (for each cycle),

- in the third column the rotation angle at which the max amplitude occurs.

this is the first lines of your file :

1.069;1.000;117.188;7.263

1.074;1.001;122.070;11.108

1.069;1.001;126.953;12.817

1.069;1.001;122.070;9.827

1.074;1.001;117.188;9.827

1.074;1.002;126.953;11.536

1.069;0.999;126.953;12.817

1.069;1.002;122.070;10.254

1.064;1.002;131.836;10.254

1.069;1.003;131.836;6.836

1.064;1.002;122.070;3.418

1.079;1.003;131.836;0.000

1.064;1.003;131.836;-8.118

1.064;1.003;131.836;-11.536

and so forth ...

maybe the numbers here are what you describe , but in that case the work is done , no ? the data are in a file and that was it...

[Michael]

Hi,

Maybe I wrote it wrong. So once again. Yes, there are 4 columns in the file, where in the fourth column is the amplitude. There are 100 cycles stored in this file one under another. Each cycle has 7200 records.

I would like to create a new matrix:

-Column 1: I want to generate the cycle number (from 1 to 100)

-Column 2: applies a low pass filter (up to 2.5 kHz) to the 4th column of the file and registers the max amplitude from 300 to 500

-Column 3: angle at which is the maximum amplitude

I hope this is now understandable.

[Mathieu NOE]

ok got it this time !

% These are data concerning engine vibrations.
% One file is one cycle. In one cycle the shaft rotates twice,
% i.e. 720 degrees, the measurement was made 7200 times per cycle,
% i.e. every 0.1 degree of rotation. 

clc
clearvars

% parameters
RPM = 2000; % my assumption
Fs = RPM/60*3600; % 3600 samples for one shaft rev
dt = 1/Fs;
nfft = 500; % frame length
overlap = 0.75; % window overlap % overlap is 95% of nfft here
angle_range = [300 500]; % range from 300 - 500 degrees, because there is the combustion 
freq_range = [0 10]; % (units : kHz !) range to display on Frequency axis (Spectrogram y axis)

% load data
out = readmatrix('4317.dat'); % 100 cycles stored in this file one under another. Each cycle has 7200 records.
[samples,cols] = size(out);
record_length = 7200;
cycles = samples/record_length;
angl_increment = 0.1;
angl = angl_increment*(0:record_length-1); % see above : One file is one cycle. In one cycle the shaft rotates twice,
                        % i.e. 720 degrees, the measurement was made 7200 times per cycle,
                        % i.e. every 0.1 degree of rotation.
                        
% low pass filter %
fc_lpf = 2500;       % LPF cut off freq
N_lpf = 2;           % filter order
[b,a] = butter(N_lpf,2*fc_lpf/Fs);
    
for ci =1:cycles
    start = 1+(ci-1)*record_length;
    stop = start + record_length -1;
    data_one_cycle = out(start:stop,:); % all 4 columns one cycle (angle from 0 to 720 by 0.1 deg increment)
    % keep only data from angle 300 to 500 deg
    ind = (angle_range(1)/angl_increment:angle_range(2)/angl_increment); 
 
    % -Column 1: I want to generate the cycle number (from 1 to 100)
    cycle_number(ci) = ci; % easy !

    % -Column 2: applies a low pass filter (up to 2.5 kHz) to the 4th column of the file and registers the max amplitude from 300 to 500 
    % -Column 3: angle at which is the maximum amplitude
    col4_lpf = filtfilt(b,a,data_one_cycle(ind,4));
    [col4_lpf_max(ci),indmax] = max(col4_lpf); % max of the positive side or max of absolute ? 
    angle_maximum_amplitude(ci) = angle_range(1) + indmax*angl_increment;
        
end

plot(cycle_number,col4_lpf_max,cycle_number,angle_maximum_amplitude);
legend('Max signal amplitude ','Angle maximum amplitude ');
xlabel('Cycle number');

data_out = [cycle_number(:) col4_lpf_max(:) angle_maximum_amplitude(:)];
% -Column 1: cycle number (from 1 to 100)
% -Column 2: applies a low pass filter (up to 2.5 kHz) to the 4th column of the file and registers the max amplitude from 300 to 500 
% -Column 3: angle at which is the maximum amplitude

[Michael]

Thank you very much for your help, the script works great!