玩陀螺儀的都會遇到一個問題就是,陀螺儀輸出的是角速度和線加速度。怎麼把加速度轉化成位移就值得研究一下。
首先我們講一下傅立葉變換,傅立葉本身就是一個線性積分變換。主要是将信号在時域和頻域中進行變換。因為我們相信任何一個信号都可以分解成sin函數。sin函數的頻率,振幅可以組合成很多的信号形式。
傅立葉變換的數學公式是這樣的。
![](https://img.laitimes.com/img/_0nNw4CM6IyYiwiM6ICdiwiIx0DciV2dmADM30zd-cGcq5CaWpXTwMGVZxmWE5UbKpmT310RapmVy0EeFdVT4NGVPpXU61EaoRVW0UFROtmSt5kaS1WWtlTeaJjTzwUeWdkW1ZVbjZ3ZHRGaxIDTox2RaxWMywEeZNDWw4EWalXOTF2dGJDTup0MiVXRXF2aWdlYwR3VhNTO5xkNNh0YwIFSh9CXt92YuM3YltWas5iclN3Ztl2Lc9CX6MHc0RHaiojIsJye.jpg)
簡單的有一個示意動畫就可以說明問題。
經過傅立葉變換,我們所謂的頻域積分也都是基于sin函數的積分。
傅立葉函數有個積分性質,當積分函數進行傅立葉變換的時候有下面這個特性
其中機關脈沖函數是這樣的
其中機關階躍函數表現是這樣的,
是以一個函數的積分的傅立葉變換應該等于
用matlab進行驗證,程式來源參考見尾部,
clc;
clear;
load(\'walk3.1.txt\');
time =[];
for i=0:1193
time = [time;i];
end
signal = sin(0.01*time*pi);
figure
plot(signal)
%y = walk3_1(:,6);
y = signal;
velocity =[];
for i=1:1194
velocity =[velocity;(sum(y(1:i)))];
end
distance = [];
for i=1:1194
distance =[distance;(sum(velocity(1:i)))];
end
figure;
subplot(2,1,1)
plot(velocity)
subplot(2,1,2)
plot(distance)
title(\'time domain - velocity -distance \')
figure;
z = fft(y);
z1 = fftshift(z);
abs1 = abs(z);
abs2 = abs1(1:1194/2+1);
abs3 = abs(z1);
f = (0:1193);
subplot(3,1,1)
plot(y)
subplot(3,1,2)
plot(abs1)
subplot(3,1,3)
plot(abs3)
title(\'fft - original - fft_abs -fftshift+abs \')
L = numel(y);
T = L;
df= 1/T;
if ~mod(L,2)
f = df*(-L/2:L/2-1);
else
f = df*(-(L-1)/2:(L-1)/2);
end
w = 2*pi*f;
for i = 1:numel(z1)
z3(i) = z1(i)*(-1i)/w(i)+z1(i);
end
yvel = ifft(ifftshift(z3),\'symmetric\');
z4 = fftshift(fft(velocity));
for i = 1:numel(z4)
z5(i) = z4(i)*(-1i)/w(i)+z4(i);
end
ydis = ifft(ifftshift(z5),\'symmetric\');
iomegaOutVel = iomega(signal,1,3,2);
figure
hold on
grid on
plot(velocity)
plot(yvel)
plot(iomegaOutVel)
title(\'Frequency domain - velocity - fft+velocity \')
figure
hold on
grid on
plot(distance)
plot(ydis)
title(\'Frequency domain - distance - fft+distance \')
iomegaOutDis = iomega(velocity,1,2,1);
plot(iomegaOutDis)
其中有一個iomega函數是一個是一個自定義的函數
function dataout = iomega(datain, dt, datain_type, dataout_type)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% IOMEGA is a MATLAB script for converting displacement, velocity, or
% acceleration time-series to either displacement, velocity, or
% acceleration times-series. The script takes an array of waveform data
% (datain), transforms into the frequency-domain in order to more easily
% convert into desired output form, and then converts back into the time
% domain resulting in output (dataout) that is converted into the desired
% form.
%
% Variables:
% ----------
%
% datain = input waveform data of type datain_type
%
% dataout = output waveform data of type dataout_type
%
% dt = time increment (units of seconds per sample)
%
% 1 - Displacement
% datain_type = 2 - Velocity
% 3 - Acceleration
%
% 1 - Displacement
% dataout_type = 2 - Velocity
% 3 - Acceleration
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Make sure that datain_type and dataout_type are either 1, 2 or 3
if (datain_type < 1 || datain_type > 3)
error(\'Value for datain_type must be a 1, 2 or 3\');
elseif (dataout_type < 1 || dataout_type > 3)
error(\'Value for dataout_type must be a 1, 2 or 3\');
end
% Determine Number of points (next power of 2), frequency increment
% and Nyquist frequency
N = 2^nextpow2(max(size(datain)));
df = 1/(N*dt);
Nyq = 1/(2*dt);
% Save frequency array
iomega_array = 1i*2*pi*(-Nyq : df : Nyq-df);
iomega_exp = dataout_type - datain_type;
% Pad datain array with zeros (if needed)
size1 = size(datain,1);
size2 = size(datain,2);
if (N-size1 ~= 0 && N-size2 ~= 0)
if size1 > size2
datain = vertcat(datain,zeros(N-size1,1));
else
datain = horzcat(datain,zeros(1,N-size2));
end
end
% Transform datain into frequency domain via FFT and shift output (A)
% so that zero-frequency amplitude is in the middle of the array
% (instead of the beginning)
A = fft(datain);
A = fftshift(A);
% Convert datain of type datain_type to type dataout_type
for j = 1 : N
if iomega_array(j) ~= 0
A(j) = A(j) * (iomega_array(j) ^ iomega_exp);
else
A(j) = complex(0.0,0.0);
end
end
% Shift new frequency-amplitude array back to MATLAB format and
% transform back into the time domain via the inverse FFT.
A = ifftshift(A);
datain = ifft(A);
% Remove zeros that were added to datain in order to pad to next
% biggerst power of 2 and return dataout.
if size1 > size2
dataout = real(datain(1:size1,size2));
else
dataout = real(datain(size1,1:size2));
end
return
其中以sin函數為輸入信号,做兩次積分得到位移
其中三個積分分别給出了三個結果,具體原因有待進一步分析。有興趣的小夥伴可以看看我的程式哪裡出現了問題。
後續有待繼續研究。
參考網站:
Matlab計算頻域積分
Matlab數值積分實作