single sided fast fourier transfrom from a data file

3 ビュー (過去 30 日間)
Hussein Kokash
Hussein Kokash 2022 年 10 月 27 日
コメント済み: Mathieu NOE 2022 年 11 月 3 日
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

採用された回答

Mathieu NOE
Mathieu NOE 2022 年 10 月 27 日
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

その他の回答 (1 件)

Hussein Kokash
Hussein Kokash 2022 年 10 月 27 日
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!
  8 件のコメント
Hussein Kokash
Hussein Kokash 2022 年 11 月 3 日
I will give it a try and get back to you!
Thanks again Mathieu, you have been a great help!
Mathieu NOE
Mathieu NOE 2022 年 11 月 3 日
as always, my pleasure !

サインインしてコメントする。

カテゴリ

Help Center および File ExchangeFourier Analysis and Filtering についてさらに検索

製品


リリース

R2018a

Community Treasure Hunt

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

Start Hunting!

Translated by