時域積分與頻域積分 實現及對比
阿新 • • 發佈:2017-12-07
wal 數學公式 save 參考 end hold on 說明 media points
玩陀螺儀的都會遇到一個問題就是,陀螺儀輸出的是角速度和線加速度。怎麽把加速度轉化成位移就值得研究一下。
首先我們講一下傅立葉變換,傅立葉本身就是一個線性積分變換。主要是將信號在時域和頻域中進行變換。因為我們相信任何一個信號都可以分解成sin函數。sin函數的頻率,振幅可以組合成很多的信號形式。
傅立葉變換的數學公式是這樣的。
簡單的有一個示意動畫就可以說明問題。
經過傅立葉變換,我們所謂的頻域積分也都是基於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數值積分實現
時域積分與頻域積分 實現及對比