FFT: extract amplitude ratio when signal is not "straight"

Question

I need to do frequency analysis on an integrating process where input is torque and output is position. If the input is a sinusoid, the output becomes something like this:

The code I use to extract amplitude ratio and phase looks like this:

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                 

If I run it on perfectly straight sinusoidal signals then it works perfectly. But as soon as I add a distortion like above, both phase and amplitude values are not right. I think I need to somehow "flatten" the signal. But I'm not sure how to extract correct amplitude from it. What is the amplitude? I would say in the output it is ~45 measured from one "plateau" to the next since that's how far the thing moves. That would be ratio of ~22.5. However the result of the calculation is 196.

Maybe I'm thinking about it wrong? I want to eventually derive a transfer function from the torque input to the positional output using experimental data. Perhaps someone can show how to do that instead?

I have been thinking that what I could do is record amplitude ratio and phase and then make a bode plot and from that easily extract the transfer function. I have so far not been able to get as far as getting a bode plot out of running tests with varying input frequencies.

Answer 1

Because FFT assumes you performing frequency analysis of perfectly periodical signal (exactly one period of the signal) your fft(out) will contain very big power disturbance (see the Periodicity and Shift theorem).

I believe, your case you can avoid of FFT analysis artifacts by performing some system modification. Instead of estimating the transfer function of the system, you can estimate the transfer function of the system + filter. I.e. you have to pass the out signal of your system through the high-pass filter:

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

Then, you can perform your analysis.

The frequency response of the filter you can estimate trough the freqz function (or just H = fft([1 -1], length(out)); in my case). Then you can eliminate the filter influence in the frequency domain by inverse response applying fft_out = fft_out ./ H(:);. Also, do not forget to nullify the 0-th frequency fft_out(1) = 0; before maximum estimation.

By the way, the phase difference estimation of the different frequencies looks strange (in your code phase = angle(fft_out(idx_out)) - angle(fft_in(idx_in))). Looks like you must use idx_in (or idx_out, depends what estimation is more reliable) for bout angles.

Note: this answer is not the complete guide and some real life enhancements possible will be required.

P.S. Try applying windowing for frequency response estimation in real world applications (for example the Hamming window).

P.P.S Try ask your question at the https://dsp.stackexchange.com/

Update: In some cases, instead of neglecting filter influence, you can perform same linear input signal transformation of the input signal: in = filter([1 -1], 1, in);