БПФ: извлечь коэффициент амплитуды, когда сигнал непрямой

Мне нужно провести частотный анализ процесса интегрирования, где вход — это крутящий момент, а выход — положение. Если вход представляет собой синусоиду, вывод становится примерно таким:

введите здесь описание изображения

Код, который я использую для извлечения соотношения амплитуд и фазы, выглядит следующим образом:

freq = 40;
freq_rad = freq * 2 * pi
phase_offset_rad = 30 * pi / 180
gain = 0                                                                                                                                                                           
fs = 500;                                                                                                                                                                          
L = 100;                                                                                                                                                                           
t = (0:L-1)*(1/fs);                                                                                                                                                                

in = 2 * sin(freq * 2 * pi * t);                                                                                                                                                   

pos_in = [];                                                                                                                                                                       
vel = 0;                                                                                                                                                                           
pos = 0;                                                                                                                                                                           
for i = 1:length(t)                                                                                                                                                                
    vel = vel + in(i);
    pos = pos + vel;
    pos_in = [pos_in; pos];                                                                                                                                                        
end                                                                                                                                                                                

out = pos_in;                                                                                                                                                                      
%out = (2 + gain) * sin(freq * 2 * pi * t + phase_offset_rad);

fft_in = fft(in);                                                                                                                                                                  
fft_out = fft(out);                                                                                                                                                                
[mag_in idx_in] = max(abs(fft_in));                                                                                                                                                
[mag_out idx_out] = max(abs(fft_out));                                                                                                                                             

phase = angle(fft_out(idx_out)) - angle(fft_in(idx_in))                                                                                                                            
phase_deg = phase / (pi / 180)                                                                                                                                                     
ratio = mag_out / mag_in                 

Если я запускаю его на идеально прямых синусоидальных сигналах, он работает отлично. Но как только я добавляю искажение, как указано выше, значения фазы и амплитуды неверны. Думаю, нужно как-то "сгладить" сигнал. Но я не уверен, как извлечь из него правильную амплитуду. Какова амплитуда? Я бы сказал, что на выходе это ~ 45 измерено от одного «плато» до другого, так как это то, как далеко движется вещь. Это будет соотношение ~ 22,5. Однако результат вычисления равен 196.

Может быть, я неправильно об этом думаю? Я хочу в конечном итоге получить передаточную функцию от входа крутящего момента к позиционному выходу, используя экспериментальные данные. Возможно, кто-то может показать, как это сделать?

Я думал, что я мог бы записать отношение амплитуд и фазу, а затем построить график Боде и легко извлечь из него передаточную функцию. До сих пор мне не удавалось получить график Боде из запуска тестов с различными входными частотами.


person Martin    schedule 05.01.2018    source источник


Ответы (1)


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

Я считаю, что в вашем случае вы можете избежать артефактов анализа БПФ, выполнив некоторую модификацию системы. Вместо оценки передаточной функции системы можно оценить передаточную функцию системы + фильтр. т.е. вы должны передать выходной сигнал вашей системы через фильтр верхних частот:

out = filter([1 -1], 1, out);

Затем вы можете выполнить свой анализ.

Частотную характеристику фильтра можно оценить с помощью функции freqz (или просто H = fft([1 -1], length(out)); в моем случае). Затем вы можете устранить влияние фильтра в частотной области с помощью обратной характеристики, применяя fft_out = fft_out ./ H(:);. Также не забудьте обнулить 0-ю частоту fft_out(1) = 0; перед максимальной оценкой.

Кстати, оценка разности фаз разных частот выглядит странно (в вашем коде phase = angle(fft_out(idx_out)) - angle(fft_in(idx_in))). Похоже, вы должны использовать idx_in (или idx_out, зависит от того, какая оценка более надежна) для углов схватки.

Примечание: этот ответ не является полным руководством, и потребуются некоторые возможные улучшения в реальной жизни.

P.S. Попробуйте применить оконный режим для оценки частотной характеристики в реальных приложениях (например, окно Хэмминга). .

P.P.S. Попробуйте задать свой вопрос на странице https://dsp.stackexchange.com/.

Обновление: В некоторых случаях вместо того, чтобы пренебрегать влиянием фильтра, вы можете выполнить такое же линейное преобразование входного сигнала: in = filter([1 -1], 1, in);

person Andrei Davydov    schedule 05.01.2018
comment
Спасибо Андрей. По крайней мере, теперь мне есть с чего начать. - person Martin; 05.01.2018
comment
Значит, для fft сумма всех элементов массива должна быть равна нулю? - person Martin; 05.01.2018
comment
Нет, в общем случае сумма отлична от нуля: sum(fft(1:16)) равно 16.0000 - 0.0000i - person Andrei Davydov; 08.01.2018