Акселерометр с FFT - странный вывод

Question

После прочтения большого количества исследований и работ по теме у меня все еще есть проблемы с применением FFT для данных акселерометра. Большая часть моего кода взята из официального примера MATLAB: FFT for one dimension. После большего чтения я нашел этот вопрос: FFT and accelerometer data: why am I getting this output?, где было предложено использовать оконный режим. Поэтому после некоторого большего чтения я добавил окно hamming к моему коду.

Мои данные выглядят, что на участке:

И это код, который я использую для быстрого преобразования Фурье:

fs = 1/0.02; %0.02 comes from picking sample each 20ms 
m = size(data,1); 
w = hanning(m); 
yw = w.*data; 
n = pow2(nextpow2(yw)); 
y = fft(yw,size(n,1)); 
f = (0:size(n,1)-1)*(fs/size(n,1)); 
power = y.*conj(y)/size(n,1); 
figure 
plot(f,power) 

Проблема заключается в том, что мой сюжет из этого кода выглядит следующим образом:

Может кто-нибудь сказать мне, что не так с моим кодом? Честно говоря, я бы исключал, что это будет выглядеть лучше (примерно так: http://imgur.com/wGs43), поэтому я задаю этот вопрос.

EDIT: Мои данные можно найти здесь: https://dl.dropboxusercontent.com/u/58774274/exp.txt

Answer 1

Ваш Fs является 50 так самая высокая частота в ваших данных может быть Fs/2 = 25Hz.

Посмотрите, подходит ли этот код.

fid = fopen('1.txt','r'); 
C = textscan(fid, '%f'); 
fclose(fid); 
data = C{1}; 
fs = 50; 
m = length(data); 
nfft = 2^nextpow2(m); 
y = fft(data,nfft)/m; 
f = fs/2 * linspace(0,1,nfft/2+1); 
power = abs(y); 
subplot(211) 
plot(f,power(1:nfft/2+1)) 
t = (0 : m-1)/fs; 
s0 = .8*fs : 3.2*fs; % .8 sec to 3.2 sec 
p(s0) = .5*cos(2*pi*3.3*t(s0)+.25*pi); 
p = p + mean(data); 
subplot(212) 
plot(t,data);hold on 
plot(t,p,'r') 

Это ваши данные в частотной области.

Существует пик в 3.3 Hz.

В качестве доказательства я график синусоидальной с частотой 3.3 Hz вместе с вашими данными,

Как вы можете видеть, это полностью соответствует вашим данным.

Answer 2

Какая ваша эталонная графика показывает только величину на логарифмической шкале (в дБ). Также кажется, что смещение постоянного тока удалено. Таким образом, вы получаете что-то вроде этого, если вы адаптируете его в своем коде.

data = data-mean(data); % DC removal 
fs = 50; 
m = length(data); 
nfft = 2^nextpow2(m); 
y = fft(data,nfft)/m; 
f = fs/2 * linspace(0,1,nfft/2+1); 
power = abs(y); 
plot(f,10*log10(power(1:nfft/2+1))); % plot log magnitude 
ylim([-30 0]) %limit axis 

По меньшей мере, похоже на меня.

Если вы дополнительно взять абсолютный квадрат вы, в основном, оценивающие спектральную плотность мощности вашего сигнала. Таким образом, вы должны изменить power = abs(y); на power = abs(y.^2);, а затем, очевидно, придется адаптировать ваш ylim([-30 0]) к чему-то ниже, чем -50.

Answer 3

Ваши участки будут выглядеть лучше, если вы сначала удалите смещение по постоянному току (вычтите среднее значение всех выборок из каждой точки перед вычислением БПФ), а затем постройте только половину точек данных результатов БПФ (N/2) или менее (так как верхняя половина результата FFT является просто сопряженным зеркалом первой половины для реального ввода данных).