single sided fast fourier transfrom from a data file

Hello everyone, I am trying to use matlab to a number of data files that contains velocity values and different positions that form a line, the velocity plot form a sinosodial shape and I would like to apply FFT and then get the frequecny, and use it to find the wavelength for each file.

(I have attached a sample file, first column position, 2nd column velcoity values, 3rd column time for the specific file)

What I expect is to convert this velocity sinosodial wave into a a single sided FFT, and then find the corresponding frequency of the highest peak which will define the wavelength at that file (labmda=1/frequency)

I have tried this but it is far from what I expect, I am supposed to get an FFT with a peak, and the 1/(corresponding x value of that peak) suppose to give the wavelength of that sinosodial wave.

data=load('data_1945.txt');

% Define x and y from the uploaded files

x = data (:,1);

y = data (:,2);

% Get Wavelength from FFT

A = fftshift(abs(fft(y)));

[pks,freqs] = findpeaks(A);

[pk1,idx1] = max(pks); %biggest peak and its index

pk_max = pk1;

idx_max = idx1;

f1 = idx_max; %frequency of biggest peak

lambda = 1 / f1;

figure(1)

plot(A)

hold on

figure(2)

plot(x,y)

hold on

Answer 1

[Mathieu NOE]

hello

if your number of periods recorded is not very high and your signal is "clean" then I would prefer to mesure the time difference between succesives zero crossing points (interpolated)

this will be definitively more accurate than a fft with only few samples (frequency resolution df = Fs / samples)

data = readmatrix('data_1945.txt');
t = data(:,1);
fs = 1/mean(diff(t));
y = data(:,2);


% fft method
[fhz,fft_spectrum] = do_fft(t(:),y(:));
[amplitude1, idx] = max(fft_spectrum);
frequency = fhz(idx);
period1 = 1/frequency
amplitude1

% zero crossing method -alternative for clean periodic signals-
zct = find_zc(t,y,0);
period2 = mean(diff(zct))
amplitude2 = 0.5*(max(y) - min(y))

figure(1)

plot(t,y,'b',zct,zeros(1,numel(zct)),'*r','markersize',25)
title('target position: 1.5 radians')
xlabel('Time[s]')
ylabel('Position[radians]')
legend('signal','zc points')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function zct = find_zc(x,y,threshold)

% positive slope "zero" crossing detection, using linear interpolation
y = y - threshold;
zci = @(data) find(diff(sign(data))>0);  %define function: returns indices of +ZCs
ix=zci(y);                      %find indices of + zero crossings of x
ZeroX = @(x0,y0,x1,y1) x0 - (y0.*(x0 - x1))./(y0 - y1); % Interpolated x value for Zero-Crossing 
zct = ZeroX(x(ix),y(ix),x(ix+1),y(ix+1));

end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [freq_vector,fft_spectrum] = do_fft(time,data)

dt = mean(diff(time));
Fs = 1/dt;
nfft = length(data); % maximise freq resolution => nfft equals signal length

% window : hanning
window = hanning(nfft);
cor_coef = length(window)/sum(window);

% fft scaling
% fft_spectrum = abs(fft(data))/nfft;
fft_spectrum = abs(fft(data.*window))*2*cor_coef/nfft;

% 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

Answer 2

[Hussein Kokash]

Hello Mathieu, thanks for your input!

Actullay at some time steps, the signal looks messy, something like this:

That is why am trying to implement FFT to get an accurate outcome of the amplitude through frequency.

My idea is to get the peak of every FFT results and find its index of frequecy, then 1/frequency to find wavelength.

Any suggestions? appreciate it!

[Mathieu NOE]

hello again

do you have a record with this kind of signal ? can you share it ?

M

[Hussein Kokash]

Yes sure, I have put two files at two different time steps.

[Hussein Kokash]

The idea is, although they are somehow not clean, but they have higher wavelength value, and at later times, the wavelength drops to lower value (splitting), supposed to get something like this:

So this is Wavelength Vs. Time, dont't mid the axis range it is not true I just wanted to give an idea of what am looking for.

[Mathieu NOE]

hello again

so I suppose you have one bigger data file from which we could compute this wavelength vs time

why are you then splitting this bigger file in smaller chuncks (the four extracts you posted) ?

maybe it's better to work on the full data in once rather than splitting it in small chuncks

[Hussein Kokash]

Hello Mathieu,

It is related to time step, each file is a certain time step, but now it is better as I remove the most noisy files and keep the ones that produces resaonable output.

Another question for you, do you have the code you posted in a for-loop for multiple files? if yes can you share it please?

Thanks again, I appreciate it!

[Mathieu NOE]

hello again

you can try this code for multiple files (for loop)

path = pwd; % put current path here 
S = dir(fullfile(path,'data_*.txt')); 


for k = 1:numel(S) 
    
    folder = S(k).folder; 
    filename = S(k).name; 
    F = fullfile(S(k).folder,S(k).name); 
    data = readmatrix(F); 
    t = data(:,1); 
    fs = 1/mean(diff(t)); 
    y = data(:,2); 
    
    
    % zero crossing method -alternative for clean periodic signals- 
    zct = find_zc(t,y,0); 
    
    disp('--------------------------'); 
    disp(['Results for file : ' filename]); 
    period2 = mean(diff(zct)); 
    amplitude2 = 0.5*(max(y) - min(y)); 
    disp(['Period : ' num2str(period2)]); 
    disp(['Amplitude : ' num2str(amplitude2)]); 
    disp('--------------------------'); 
    
    figure(k) 
    
    plot(t,y,'b',zct,zeros(1,numel(zct)),'*r','markersize',25) 
    filename_title = filename; 
    filename_title = strrep(filename_title,'_',' '); % replace underscore with space tab 
    title(filename_title) 
    xlabel('Time[s]') 
    ylabel('Position[radians]') 
    legend('signal','zc points') 
    
end 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function zct = find_zc(x,y,threshold) 

% positive slope "zero" crossing detection, using linear interpolation 
y = y - threshold; 
zci = @(data) find(diff(sign(data))>0); %define function: returns indices of +ZCs 
ix=zci(y); %find indices of + zero crossings of x 
ZeroX = @(x0,y0,x1,y1) x0 - (y0.*(x0 - x1))./(y0 - y1); % Interpolated x value for Zero-Crossing 
zct = ZeroX(x(ix),y(ix),x(ix+1),y(ix+1)); 

end 

[Hussein Kokash]

I will give it a try and get back to you!

Thanks again Mathieu, you have been a great help!

[Mathieu NOE]

as always, my pleasure !