필터 지우기
필터 지우기

single sided fast fourier transfrom from a data file

조회 수: 2 (최근 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 CenterFile Exchange에서 Fourier 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