Discrete input signal and continuous transfer function

I have a discrete input signal, something like an earthquake ground motion, measured every few milliseconds, with a specified units. That input signal has been converted to frequency-domain through FFT.

I also have a trasnfer function which I need to run the discrete input signal through while in frequency domain (I assume), which is very complex. I have succeffully modeled and plotted it with the tf function in MATLAB.

Is there a way to multiply the discrete input signal with the transfer function? I have tried the function lsim function, but it said that a continous system is required. Should I covert the transfer function to a discrete function and the multiply it?

Thank You so much for the help

[Sam Chak]

@Nerma Caluk, do you mean that you want to inject the Amplitude Spectrum (in Freq domain) of discrete-time input signal S into a continuous-time transfer function?

Fs = 1000;            % Sampling frequency                    
T = 1/Fs;             % Sampling period       
L = 1500;             % Length of signal
t = (0:L-1)*T;        % Time vector

S = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);
stem(S(1:50))
title("Discrete-time input signal, S(n)")
xlabel("n")
ylabel("S(n)")

Y = fft(S);
P2 = abs(Y/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);

f = Fs*(0:(L/2))/L;
plot(f,P1) 
title("Single-Sided Amplitude Spectrum of S(t)")
xlabel("f (Hz)")
ylabel("|P1(f)|")

[Nerma Caluk]

Yes, that is exaclty what I meant! I have managed to convert the trasnfer function to discrete system with c2d function, and then use the lsim function. The lsim function did not require amplitude spectrum of the signal, just the original signal and discretized transfer function, but it is giving me very high units. I am not completley sure what exaclty happens in the lsim function and I see no useful examples online.

Answer 1

[Star Strider]

It would likely be best to convert the continuous-time transfer function to a discrete-time transfer function (ideally using the Tustin transformation, using the sampling frequency of the earthquake ground motion signal as the sampling frequency of the discrete filter) and then use that to filter the signal. If you have the Signal Processing Toolbox, this should be relatively straightforward. Realise the discrete transfer function as a second-order-section realisation (rather than a transfer function realisation) for the best results.

[Nerma Caluk]

Thank You for helping. This is part of the code that I have developed so far, excluding the exctration of the data part since it is a bit lengthy. The code contains all explanation, but I do not know now how to divide/multiply (not sure which at the moment) the Amplitude Spectra after FFT on SPZ with the Transfer function which I have convreted into discrete function thorugh function c2d and used Tustin transformation. Any idea how to proceed?

clear

clc

close all

% Extraction of the data was detelted from this code since it is very long

% The extracte signal data is identified as SPZ which is sampled every

% 0.01887 sec

% SPZ - discrete values of change of signal (specific units)

% ms - correspoindng time values in milliseconds from the input data

Fs = 53; %Sampling Frequency (Hz)

T = 1/Fs; %Sampling Period (sec)

L = length(ms); %Length of signal

t = ms./1000; %Time vector

omega_NyqS = pi/T;

freq_NyqS = omega_NyqS/(2*pi);

% plotting of original signal

figure

plot(t,SPZ,'Line Width', 3);

title('Original signal before filtering')

xlabel('t (Seconds)')

ylabel('Digital Units')

% Computation of FFT of the input signal

Y = fft(SPZ);

Y = Y./L;

figure

plotPowerSpectrum(Y, Fs,'Original');

% Half of the Y array is positive frequencies

f = Fs*(0:length(Y)/2)/(length(Y)/2);

w = 2*pi*f;

% Defining the transfer function

R_g = 1800; % coil resistance [ohms]

R_s = 2680; % damping resistance [ohms]

G_1 = 175; % generator constant of the magnet-coil system [V/(m/s^2)]

G_2 = 23700; % pre-amplifier gain

omega_b = 0.31416; % output high pass cut-off angular frequency (rad./s)

omega_p = 57.1199; % output low pass cut-off angular frequency (rad./s)

f_s0 = 1; % responant frequency of the pendulum (Hz)

h = 0.85; % damping constant

S_R = 53; % sampling rate, nominal (Hz)

omega_0 = 2*pi*f_s0;

K = 204.8;

s = tf('s');

G = R_s/(R_g + R_s);

S_p = s/(s^2 + 2*h*omega_0*s + omega_0^2);

num = omega_p.^2;

denom_1 = s^2 + 2*cos(pi/8)*(omega_p)*s + (omega_p)^2;

denom_2 = s^2 + 2*cos((3*pi)/8)*(omega_p)*s + (omega_p)^2;

F_sa = num/denom_1;

F_sb = num/denom_2;

F_s = [(F_sa)^2]*[(F_sb)^2];

F_b = s/(omega_b + s);

A_SP_DU = K*G*G_1*G_2*S_p*F_b*F_s; %Units [V/(m/s^2)]

Tf = A_SP_DU;

Tfz = c2d(Tf,0.001,'tustin');

% Theoretically, now I need to divide the Amplitude Spectra (FFT values of

% Y) with the Transfer Fucntion?

function plotPowerSpectrum(y,fs,name)

n = length(y); % number of samples

f = (0:n-1)*(fs/n); % frequency range

power = abs(y).* abs(y); % power of the DFT

loglog(f,power)

title(name);

xlabel('Frequency')

ylabel('Power')

grid

drawnow

end

Thank You all for such quick help and advices!

[Sam Chak]

Hi @Nerma Caluk

ms is missing. Any idea what the value is?

% ms - correspoindng time values in milliseconds from the input data
Fs = 53;        %Sampling Frequency (Hz)
T = 1/Fs;       %Sampling Period (sec)
L = length(ms); %Length of signal 

Unrecognized function or variable 'ms'.

t = ms./1000;   %Time vector

[Nerma Caluk]

SPZ and ms are the only 2 input values from the gathered data. ms is supposed to be a time vector with units of milliseconds, gatherede every 0.0188679 seconds. SPZ is the y-axis data that shows how the signal changes.

Their extraction is bit lengthy code and lead to array of 300005 data points, and I took it out but any arbitay values should work.

[Star Strider]

I would just use lsim to do the filtering, however the Bode plot of the filter seems that it does not do much actual filtering —

% Extraction of the data was detelted from this code since it is very long
% The extracte signal data is identified as SPZ which is sampled every
% 0.01887 sec
% SPZ - discrete values of change of signal (specific units)
% ms - correspoindng time values in milliseconds from the input data

SPZ = randn(1, 53*100);                                   % Create Data
ms = linspace(0, numel(SPZ)-1, numel(SPZ))/53;     % Create Data
SPZ = SPZ + 10*sum(sin([1;150;200]*2*pi*ms/numel(ms)));

figure
plot(ms, SPZ)

Fs = 53;        %Sampling Frequency (Hz)
T = 1/Fs;       %Sampling Period (sec)
L = length(ms); %Length of signal 
t = ms./1000;   %Time vector
omega_NyqS = pi/T;
freq_NyqS = omega_NyqS/(2*pi);

% plotting of original signal
figure
plot(t,SPZ,'LineWidth', 3);
title('Original signal before filtering')
xlabel('t (Seconds)')
ylabel('Digital Units')

% Computation of FFT of the input signal
Y = fft(SPZ);
Y = Y./L;

figure
plotPowerSpectrum(Y, Fs,'Original');

% Half of the Y array is positive frequencies
f = Fs*(0:length(Y)/2)/(length(Y)/2);
w = 2*pi*f;
% Defining the transfer function
R_g = 1800;             % coil resistance [ohms]
R_s = 2680;             % damping resistance [ohms]
G_1 = 175;              % generator constant of the magnet-coil system [V/(m/s^2)]
G_2 = 23700;            % pre-amplifier gain
omega_b =   0.31416;    % output high pass cut-off angular frequency (rad./s)
omega_p =   57.1199;    % output low pass cut-off angular frequency (rad./s)
f_s0 =       1;          % responant frequency of the pendulum (Hz)
h =         0.85;       % damping constant
S_R = 53;               % sampling rate, nominal (Hz)
omega_0 = 2*pi*f_s0;
K = 204.8;
s = tf('s');
G = R_s/(R_g + R_s);
S_p = s/(s^2 + 2*h*omega_0*s + omega_0^2);
num = omega_p.^2;
denom_1 = s^2 + 2*cos(pi/8)*(omega_p)*s + (omega_p)^2;
denom_2 = s^2 + 2*cos((3*pi)/8)*(omega_p)*s + (omega_p)^2;
F_sa = num/denom_1;
F_sb = num/denom_2;
F_s = [(F_sa)^2]*[(F_sb)^2];
F_b = s/(omega_b + s);
A_SP_DU = K*G*G_1*G_2*S_p*F_b*F_s;  %Units [V/(m/s^2)]
Tf = A_SP_DU;
Tfz = c2d(Tf,T,'tustin')

Tfz =
 
  2031 z^11 + 1.421e04 z^10 + 3.858e04 z^9 + 4.264e04 z^8 - 1.218e04 z^7 - 8.529e04 z^6 - 8.529e04 z^5 - 1.218e04 z^4 + 4.264e04 z^3 + 3.858e04 z^2 + 1.421e04 z + 2031
  ---------------------------------------------------------------------------------------------------------------------------------------------------------------------
                z^11 - 5.707 z^10 + 15.19 z^9 - 25.07 z^8 + 28.53 z^7 - 23.48 z^6 + 14.22 z^5 - 6.316 z^4 + 2.009 z^3 - 0.4348 z^2 + 0.05772 z - 0.00359
 
Sample time: 0.018868 seconds
Discrete-time transfer function.

figure
bodeplot(Tf)

figure
bodeplot(Tfz)

% Theoretically, now I need to divide the Amplitude Spectra (FFT values of
% Y) with the Transfer Fucntion?

SPZ_Filt = lsim(Tfz, SPZ, ms);

figure
plot(ms,SPZ_Filt)
grid

function plotPowerSpectrum(y,fs,name)
n = length(y);          % number of samples
f = (0:n-1)*(fs/n);     % frequency range
power = abs(y).* abs(y);    % power of the DFT
loglog(f,power)
title(name);
xlabel('Frequency')
ylabel('Power')
grid
drawnow
end

I am not certain what to make of that filter, whether implemented as an analog or digital filter.

.

[Nerma Caluk]

So I do not need the FFT of the input signal at all?

The thing is that the input of trasnfer function is indicated to be s(omega), which is equivalent j*omega, and I have assumed I have do everything in frequnecy domain. Based on your recommendations I have to use the time-domain input values (SPZ) into the lsim function, not Y?

I apologize if I made wrong assumptions, I am learning as I go.

Thank You!

[Star Strider]

‘So I do not need the FFT of the input signal at all?’

No. Not unless you want to use it to design the filter, an approach I generally recommend. It is not necessary to use it as part of the signal processing procedure otherwise.

‘Based on your recommendations I have to use the time-domain input values (SPZ) into the lsim function, not Y?’

That is what I would do. I got lsim to work with your filter, however the filter does not do much filtering, as I noted. (The purpose of converting it to a discrete fitler was to use it with Signal Processing Toolbox filter functions.) It is possible to use it with the Control System Toolbox lsim function without doing the conversion, as was done here, since lsim talkes care of all those details itself.

.

[Nerma Caluk]

Thank You very much for the advices and help. The lsim works on the code and gives me reasonable plots, however the units of y-axis seem to be weird. The trasnfer function is specifically designed for the signals I am running through and it should make sense once all the data has been inputted.

The input signal data SPZ has units know as "Digital Units" aka DU, and the final transfer function A_SP_DU has given units [DU/(m/s^2)]. My assumption was that I need to divide the input values by the transfer function to conduct deconvolution, but I wasnt sure.

I will have to talk to one of my superiors about this process to understand why was this happening.

Thank You again!

[Star Strider]

My pleasure!

[Nerma Caluk]

Hi! I have one more question since I cannot find answers anywhere online. If my discrete original signal has certain units in y-axis, let’s say “Digital Units”, and the my discretized Transfer Function has another units, let’s say Digital Units per m/s^2, what would be the output values when the function lsim is used? Would it be the same as the input discretized signal (Digital Units), or the units of transfer functions, or something else completely? Thank You again so much!

Answer 2

[Paul]

Not sure what the issue was with lsim without seeing the code.

Is the transfer function contiuous time, e.g. a model of an analog device or physical structure excited by the earthquake? If so, lsim seems like it could be the way to go for this problem, taking advantage of the foh method if linear interpolation between sampled measurements is a good model of earthquake ground motion (assuming that the dynamics represented by the tf are slow relative to the measurement sample period). If not, more information on what the tf represents would be helpful.

For example using lsim ...

Generate some data at sample period of 1 ms

t = 0:.001:2;
equake = sin(2*pi/3*t);

Continous time transfer function

H = tf(100,[1 2*.7*10 100]);

The output

y = lsim(H,equake,t,'foh');  % works fine
plot(t,equake,t,y)

[Nerma Caluk]

The signal I am running through the trasnfer function has been convreted to frequency domain using FFT, and lsim is giving me an error:

"In time response commands, the time vector must be real, finite, and must contain monotonically increasing and evenly spaced time samples."

I have shared my code in one of the answers/comments above!

Thank You!