How to load a specific column from files + FFT function

2 views (last 30 days)
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

Accepted Answer

Mathieu NOE
Mathieu NOE on 10 Jun 2022
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

More Answers (2)

Mathieu NOE
Mathieu NOE on 7 Jun 2022
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 !
  13 Comments
Michael
Michael on 15 Jun 2022
Thank you once again for your help. I will try to do the rest of the project myself:)

Sign in to comment.


Michael
Michael on 20 Jun 2022
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?
  8 Comments
Mathieu NOE
Mathieu NOE on 24 Jun 2022
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
Michael on 27 Jun 2022
Thank you very much for your help, the script works great!

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!