Thank you for your answer!

# How to load a specific column from files + FFT function

92 views (last 30 days)

Show older comments

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

##### 0 Comments

### More Answers (2)

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

### See Also

### Community Treasure Hunt

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

Start Hunting!