积分过程频域分析：非正弦输出下FFT幅相提取问题与传递函数推导

Hey Martin, let's work through your problem with the torque-to-position integration frequency analysis and transfer function derivation. I see a couple of key issues in your code and approach that are throwing off your results, plus a clear path forward to fix things and get your Bode plot and transfer function sorted.

Your current approach uses max(abs(fft_out)) to find the dominant frequency component, but when your output distorts (from nonlinearities or incorrect integration), the signal gains harmonic components (integer multiples of the input frequency). The FFT's maximum amplitude might come from one of these harmonics instead of the fundamental input frequency—this is why your magnitude ratio and phase calculations are wrong.

Instead of hunting for the maximum amplitude, directly target the FFT bin corresponding to your input frequency. Here's how to calculate that index (MATLAB uses 1-based indexing for FFT results):

• Frequency resolution of your FFT: df = fs / L
• Index for your target frequency: target_bin = round(freq / df) + 1

Then extract the magnitude and phase only from this bin. This ignores harmonic distortion and focuses on the linear system response we care about.

Looking at your code, your integration steps don't account for the sampling time 1/fs:

vel = vel + in(i);  % Wrong: this assumes Δt=1, not 1/fs
pos = pos + vel;

Discrete integration requires multiplying each sample by the time step Δt = 1/fs to approximate the area under the curve. Without this, your position output is massively scaled up (which explains your 196 ratio instead of ~22.5). Here's the corrected integration:

dt = 1/fs;
vel = vel + in(i)*dt;
pos = pos + vel*dt;

Here's your code with both fixes applied, plus clearer handling of FFT components:

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;
dt = 1/fs;  % Add time step for correct integration
for i = 1:length(t)
    vel = vel + in(i)*dt;
    pos = pos + vel*dt;
    pos_in = [pos_in; pos];
end
out = pos_in;

% Calculate FFT and target fundamental frequency bin
fft_in = fft(in)/L;  % Normalize FFT by number of samples
fft_out = fft(out)/L;
df = fs/L;
target_bin = round(freq / df) + 1;

% Extract magnitude and phase for fundamental frequency
mag_in = 2*abs(fft_in(target_bin));  % *2 for single-sided amplitude
mag_out = 2*abs(fft_out(target_bin));
phase_rad = angle(fft_out(target_bin)) - angle(fft_in(target_bin));
phase_deg = phase_rad * (180/pi);
ratio = mag_out / mag_in;

% Print results
disp(['Magnitude Ratio: ', num2str(ratio)]);
disp(['Phase Shift (deg): ', num2str(phase_deg)]);

Note: We normalize the FFT by L and multiply by 2 to get the actual time-domain amplitude (since MATLAB's FFT returns a double-sided spectrum).

To get your torque-to-position transfer function, follow these steps:

• Run frequency sweep tests: For a range of frequencies (e.g., 1Hz to 100Hz), input a steady-state sinusoidal torque, wait for transient effects to die down, then collect input and output data.
• Extract fundamental components for each frequency: Use the corrected FFT method above to get magnitude ratio and phase shift for each test frequency.
• Build the Bode plot: Convert magnitude ratios to decibels (20*log10(ratio)) and plot against log-frequency; plot phase shift (in degrees) against log-frequency.
• Fit the transfer function: In MATLAB, use system identification tools:
  1. Package your data into an iddata object:

     data = iddata(out, in, dt);
     

  2. Fit a transfer function (your system is theoretically a double integrator, so try a 2nd-order model):

     sys = tfest(data, 2);  % 2nd-order transfer function
     

  3. Compare the fitted system's Bode plot with your experimental data to validate.

If your output has distortion (e.g., from mechanical saturation or friction), focusing on the fundamental frequency component (as we did) is still valid for extracting the linear part of your system's transfer function. Nonlinearities create harmonics, but they don't affect the linear system's response at the input frequency. For extra robustness, you could add a narrow band-pass filter around the input frequency before processing the output signal.

内容的提问来源于stack exchange，提问作者Martin