Why do I need to multiply the frequencies with a number, to get correct "shift" in the bode plot?

Question

Assume that we got a sine wave function

$$u(t) = A\sin(2\pi \omega(t)t)$$

Where the frequency $\omega(t)$ changes over time $t$ and $A$ is the amplitude.

Assume that we apply that $u(t)$ signal onto a transfer function $G(s)$ and we receive the output $y(t)$ in the time domain.

Then our goal is to create the bode plot diagram of $u(t)$ and $y(t)$.

I begin to do the FFT analysis.

$$F_u(t) = \mathcal{F.T}\bigl\{u(t)\bigr\}$$ $$F_y(t) = \mathcal{F.T}\bigl\{y(t)\bigr\}$$

Then I compute the ratio

$$H(t) = \frac{F_y(t)}{F_u(t)}$$

And then I cut everything in half

H = H(1:end/2)
w = w(1:end/2)

And then I plot the diagram where the y-axis is $20\log_{10}(| H |)$ and x-axis is w

Problem:

The problem is that the frequencies are not correct. To "shift" them to the correct values , I need to multiply with a constant a:

w = w*a

Question:

How do I find that constant a?

How to reproduce the issue:

close all
% Input and model
N = 30000;
t = linspace(0.0, 50, N);
w = linspace(0, 100, N);
u = 10*sin(2*pi*w.*t);
G = tf([3], [1 0.5 30]);

% Do frequency response
y = lsim(G, u, t);
close all

% Do FFT
fy = fft(y);
fu = fft(u);
H = fy./fu;

% Windowing - Half
H = H(1:end/2);
a = 1;
w = w(1:end/2)*a;

semilogx(w, 20*log10(abs(H))); 

% Do bode without phase
bode(G);

Goal:

My goal/aim is to find the transfer function from measurement data. Estimation/identification.

Answer 1

There are a couple of issues here.

1. Have a look at how to properly generate a linear chirp (https://en.wikipedia.org/wiki/Chirp)
2. Make sure to choose a frequency range consistent with your sampling frequency (i.e. from $0$ to $f_S/2$)
3. Generate the frequency axis for the plot taking into account the FFT size and the sampling frequency (also note that bode uses the angular frequency in rad/s for the x-axis, so you need to multiply the frequency in Hertz by $2\pi$)

close all
% Input and model
N = 30000;
t = linspace(0.0, 50, N).';
T = t(end)-t(1); % chirp duration

fs = 1/(T/N);  % sampling frequency
f0 = 0;        % start frequency
f1 = fs/2;     % final frequency
c = (f1-f0)/T; % chirp rate

u = 10*sin(2*pi*(c/2*t.^2 + f0*t));

G = tf([3], [1 0.5 30]);

% Do frequency response
y = lsim(G, u, t);

Nfft = 32768; % taking the next power of 2

% Do FFT
fy = fft(y,Nfft);
fu = fft(u,Nfft);
H = fy./fu;

% (Angular) frequency axis for the plots
omega = [0:Nfft-1]/Nfft*2*pi*fs;

% Windowing - Half
H = H(1:end/2);
omega = omega(1:end/2);

% Do bode without phase
bode(G);

% Overlay the curves on the same plots
subplot(211); 
hold on;
semilogx(omega, 20*log10(abs(H)),'r--'); 
subplot(212); 
hold on;
semilogx(omega, (angle(H)/pi*180),'r--'); 

Answer 2

Note this part of your code:

t = linspace(0.0, 50, N);
w = linspace(0, 100, N);
u = 10*sin(2*pi*w.*t);

The instantaneous frequency at time index $n$ is given by the phase difference between points $n-1$ and $n$.

In your case, the phase difference is $2\pi \Delta (wt)$, and this is something like $2\pi(w\Delta t + t\Delta w + \Delta t \Delta w)$, which is not what you want. You want something like $w \Delta t$.

In order to make proper frequency modulation, you have to integrate instantaneous frequency in order to get the proper phase at each time index.

Something like (in Octave code):

# accumulator
bac = 1; aac = [1 -1];
phases = filter(bac, aac, w);
x = sin( phases );

where the vector $w$ contains the instantaneous frequency (in radians per sample).

The above needs correction for units. Your code does not properly show sampling frequency, so you must transform radians and samples to your time and frequency units.