How to load a specific column from files + FFT function

2 次查看（过去 30 天）
Michael
Michael 2022-6-7
评论： Michael 2022-6-27
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

请先登录，再回答此问题。

采纳的回答

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

更多回答（2 个）

Mathieu NOE
Mathieu NOE 2022-6-7
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 个评论
显示 11更早的评论
Michael
Michael 2022-6-15
Thank you once again for your help. I will try to do the rest of the project myself:)

请先登录，再进行评论。

Michael
Michael 2022-6-20
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 个评论
显示 6更早的评论
Mathieu NOE
Mathieu NOE 2022-6-24
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 2022-6-27
Thank you very much for your help, the script works great!

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Transforms 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by