Discrete input signal and continuous transfer function

Question

Nerma Caluk 2022-10-14

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1826763-discrete-input-signal-and-continuous-transfer-function

评论： Nerma Caluk 2022-10-18

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

2 个评论
显示无隐藏无

Sam Chak 2022-10-15

在 MATLAB Online 中打开

@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 2022-10-18

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 2022-10-15

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1826763-discrete-input-signal-and-continuous-transfer-function#answer_1075298

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.

9 个评论
显示 7更早的评论隐藏 7更早的评论

Nerma Caluk 2022-10-15

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!

Star Strider 2022-10-15

在 MATLAB Online 中打开

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.

.

Star Strider 2022-10-15

My pleasure!

Nerma Caluk 2022-10-18

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 2022-10-15

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1826763-discrete-input-signal-and-continuous-transfer-function#answer_1075973

在 MATLAB Online 中打开

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)

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Nerma Caluk 2022-10-15

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!

请先登录，再进行评论。

Discrete input signal and continuous transfer function

2 个评论
显示无隐藏无

采纳的回答

9 个评论
显示 7更早的评论隐藏 7更早的评论

更多回答（1 个）

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

Discrete input signal and continuous transfer function

2 个评论 显示 无隐藏 无

采纳的回答

9 个评论 显示 7更早的评论隐藏 7更早的评论

更多回答（1 个）

1 个评论 显示 -1更早的评论隐藏 -1更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

2 个评论
显示无隐藏无

9 个评论
显示 7更早的评论隐藏 7更早的评论

1 个评论
显示 -1更早的评论隐藏 -1更早的评论