# adding subtracting PSD vibration data

We performed some vibration tests on a shaker, with two accelerometers on our setup. One was on the shaker table, to measure the input signal. And the other one was on a flexible part of our test setup, this is considered the output signal. The input was zero mean band limited white noise acceleration, between 100 and 1200 Hz.

The data that we have from the accelerometers are the acceleration spectral densities, in g^2/Hz, of the input and output. These are the absolute accelerations of the input and output. But now I need the relative acceleration, between output and input. And also the relative velocity and the relative displacement between input and output.

Can I just add or subtract PSD values to find the relative signals?

• you mean you can divide them, not subtract those linear quantities, right? (I'm assuming your system is linear, i.e. if you scale the input by a factor of 2, the output will be scaled with exactly 2 as well, right?) Mar 26, 2023 at 11:47
• Unfortunately the system is not linear. That is what I am studying now. So we repeated the test at different input levels, 2 g_rms, 4 g_rms, 6 g_rms etc. Then we noticed that the output levels were nonlinear, they did not scale with the input levels. Mar 26, 2023 at 12:31
• We can write the (mechanical) stiffness and (mechanical) damping as a function of the input accelerations that I mentioned above, but it is required to plot stiffness against the displacement occurring in the system and to plot damping against velocity occurring in the system. For example, a spring force normally is plotted against the displacement difference between both ends of the spring. That requires a subtraction of the displacement of both ends if both ends are moving. Mar 26, 2023 at 12:31
• Ah well, but if its not linear, then the absolute difference isn't much help either. You will first need to select some nonlinear model of your device and try to fill in the parameters of that! Mar 26, 2023 at 12:32
• The output data looks somewhat like the response of a single mass-spring-damper (MSD) system. So I could fit the output data to the response of a MSD system. Then use that to determine the rms velocity and rms displacement. This should be done for each test level separately, given the nonlinearity. But I figured it would be more neat to get data like the rms velocity and rms displacement straight out of the measured data instead of fitting the acceleration results to a function and then use that function for estimating the rms velocity and rms displacement. Mar 26, 2023 at 12:54

The way I understand the question:

1. You have the power spectral density estimates $$S_\ddot{x}$$ and $$S_\ddot{u}$$.
2. You want to obtain the power spectral density estimate $$S_z$$ where $$z = x - u$$.
3. You would like to find $$S_x$$ and $$S_u$$ in the frequency domain.

The power spectral density of $$z$$ with $$z = x - u$$ would be the Fourier transform of the autocorrelation $$R_z(\tau)$$, which is something like (in the stationary case) $$R_z(\tau) = E\left[ \left( x(t) - z(t) \right) \left( x(t + \tau) - z(t + \tau) \right) \right] = R_x(\tau) - R_{xu}(\tau) - R_{ux}(\tau) + R_u(\tau)$$ implying $$S_z(j \omega) = S_x(j \omega) - 2 \mathrm{Re}\left( S_{xu}(j \omega) \right) + S_u(j \omega)$$ If $$x$$ and $$u$$ are uncorrelated, then $$S_z(j \omega) = S_x(j \omega) + S_u(j \omega)$$ but that is not the case for you; you need the cross spectral density function as well.

Time-integration in the frequency domain is equivalent to multiplying with $$1/j\omega$$. As it appears you are only working with power spectral density estimates, I am a bit at loss for finding good numerical methods, but I've included a MATLAB example of a crude method below:

clear

Fs = 1e5; % sampling freq.
Ts = 1/Fs;
Nsamp = 1e5; % number of samples
t = (0:Nsamp-1)*Ts; % time vector

x = randn(1,Nsamp); % excitation signal

[B,A] = butter(3,2*pi*1e4,'s'); % example dynamical system
W = tf(B,A); % plant/dynamics/filter
Wd = W*tf('s')^2; % 2nd derivative of filter output (should be at least proper for numerical reasons)

u = lsim(W,x,t,'foh'); % find the filter response (first order hold for better interpolation)
u = u(501:end); % remove most of the transient
ud = lsim(Wd,x,t,'foh'); % find the filter response
ud = ud(501:end); % remove most of the transient

WIN = kaiser(numel(u)/10,38); % windowing to reduce artifacts
[Pu,Fu] = pwelch(u,WIN,[],numel(u),Fs,'power','onesided'); % PSD of filter
[Pud,Fud] = pwelch(ud,WIN,[],numel(ud),Fs,'power','onesided'); % PSD of deriv. filter

figure(1), clf
semilogy(Fu,Pu)
hold on
semilogy(Fud,Pud)

II = (1/(2*pi))./(1i.*Fud(2:end)); % 1/(j 2 pi f)
Pudi = (II.*conj(II)).^2.*Pud(2:end); % basic time integration in freq. domain
semilogy(Fud(2:end),Pudi)


In my experience, given that you know the excitation force $$f$$ and output $$x$$, it is usually the transfer-function $$x(j\omega)/f(j\omega) = G(j \omega)$$ that is interesting. This transfer-function can be estimated using e.g. $$G(j\omega) = \frac{S_{xx}(j\omega)}{S_{xf}(j\omega)}$$ A numerical example is given below in MATLAB:

clear

Fs = 1e4; % sampling freq.
Ts = 1/Fs;
Nsamp = 1e6; % number of samples
t = (0:Nsamp-1)*Ts; % time vector

f = randn(1,Nsamp); % excitation signal

G = rss(3); % random plant/dynamics/filter

x = lsim(G,f,t,'zoh'); % find the plant response
x = x(1001:end); % remove most of the transient
f = f(1001:end);

WIN = hamming(numel(x)/4);
[Pxx,Fxx] = pwelch(x,WIN,[],numel(x),Fs,'onesided'); % Power Spectral Density
[Pxf,Fxf] = cpsd(f,x,WIN,[],numel(x),Fs,'onesided'); % Cross Power Spectral Density

Gest = Pxx./Pxf; % estimated response
[MAG,PHASE] = bode(G,2*pi*Fxx); % "true" response

%%
figure(1), clf

subplot(2,1,1)
loglog(Fxx,abs(Gest))
hold on
loglog(Fxx,MAG(:))

subplot(2,1,2)
semilogx(Fxx,(180/pi)*unwrap(angle(Gest)))
hold on
semilogx(Fxx(:),PHASE(:))