3
$\begingroup$

I am a chemical engineer so I'm not very familiar with the methods of signal processing. Nevertheless I need to apply them now on my thesis work. Most recently a question occurred to me, that is probably answered easily with a certain insight on signal processing and FFT.

Here is my problem: I am researching residence time distributions ($\text{RTD}$) on a reactor. Therefore I inject a tracer substance and measure its concentration course over time on the inlet and outlet. So I have two signals $C_{in}$ and $C_{out}$, which are related to each other by $$C_{out} = \text{RTD} * C_{in}$$ where $*$ stands for convolution.

To calculate the desired residence time distribution of the run one has to deconvolve the signal $C_{in}$ from the $C_{out}$ signal.

To achieve that I convert the two signals into the frequency domain via FFT: \begin{align} &C_{out}(t) \xrightarrow{FFT} C_{OUT}(f)\\ &C_{in}(t) \xrightarrow{FFT} C_{IN}(f) \end{align}

After that I divide $C_{OUT}$ by $C_{IN}$ since a convolution in the time domain is a multiplication in the frequency domain. After inverse FFT of the output I get the $\text{RTD}$.

Note: Because the amount of tracer that enters the reactor is the same that exits the reactor, the integrals of $C_{in}$ and $C_{out}$ with respect to $dt$ are of same value. Therefore the integral of the $\text{RTD}$-curve with respect to $dt$ has to be $1$.

Because the corresponding time values of the measurements of $C_{in}$ and $C_{out}$ aren't uniformly spaced I resample both with the sample frequency $f_s$ to get uni-distant time steps.

Now to my problem/question: I have noticed that the area under the RTD-curve that should be always $1$ is dependent on the sample frequency by: $$\text{area} = \frac{1}{f_s}$$ So if the sample times are $0.2\text{s}$ apart, the $\text{area} = 0.2$. If they're $2\text{s}$ apart, $\text{area} = 2$. Only with $1\text{s}$ steps I get $\text{area} = 1$. Why?

Thank you in advance for your answers. Here is also a section of my Matlab code:

% resampling fs=1 [Cout, t]=resample(Cout_orig, time_orig, fs); Cin = resample(Cin_orig, time_orig, fs); % deconvolution COUT=fft(Cout); CIN=fft(Cin); RTD=COUT/CIN; rtd=ifft(RTD); rtd=rtd(:,1); % calculating area area=trapz(t, rtd);

Integrals of original signals (_orig), resampled signals (_resample), fft-converted & back-converted signals (_ifft) divided by each other (see legend) and Integral of RTD over sample frequency

Improve this question
$\endgroup$
5
  • $\begingroup$ I'd start by verifying that the area under the original signal and resampled signal is the same when using trapz. If it is, there might be something going on with the FFT scaling... $\endgroup$ Commented Nov 3, 2022 at 17:14
  • $\begingroup$ As a side note, unless your input and output signals are completely noise-less, I wouldn't simply divide both FFTs to get the frequency response. With noisy measurements, you can get a better estimate using Power Spectral Densities. $H(f) = C_{yx} / P_{xx}$ where $C_{xy}$ is the Cross-Power Spectral Density and $P_{xx}$ is the Power Spectral Density (input $x$, output $y$). Matlab has functions to compute each: --continue-- $\endgroup$ Commented Nov 3, 2022 at 17:48
  • $\begingroup$ --continued-- cpsd, pwelch, or even one that does the work for you: tfestimate $\endgroup$ Commented Nov 3, 2022 at 17:48
  • $\begingroup$ @Jdip thank you for your comment. I checked it, they're the same as you can see from the added graph. The problem obviously occures at the division RTD=COUT/CIN. $\endgroup$ Commented Nov 3, 2022 at 20:41
  • $\begingroup$ is there any way for you to share .csv of your signals? $\endgroup$ Commented Nov 4, 2022 at 13:19

2 Answers 2

Reset to default
2
$\begingroup$

Therefore the integral of the RTD-curve with respect to dt has to be 1.

Are you sure of this? I understand that this comes from the theory of continuous-time systems. But you have effectively converted your continous-time system into a a discrete-time one. Based on the definition of discrete convolution

$ y[n] = \sum\limits_{k=-\infty}^{+\infty} x[k]h[n-k], $

I would rather say that the sum of the impulse response taps $\sum_n h[n]$ has to be 1.

Also, I see two potential problems in your MATLAB code:

  1. RTD=COUT/CIN; is wrong. If you want to perform element-wise division, you need RTD=COUT./CIN; (Notice the point before the division operator).
  2. Why do you need rtd=rtd(:,1);? In the best case (if rtd is a column) it will be unnecessary. In the worst case (if rtd is a row), it will only retain the first sample of rtd.

EDIT

I think I got a somewhat acceptable result with your second signal set. I set the sampling interval to the average sampling interval of your irregularly spaced samples, which results in 0.3113 seconds.

close all; clc; clear; data = readmatrix("Signals2.csv"); % Read skipping the first sample t_raw = data(2:end,1); x_raw = data(2:end,2); y_raw = data(2:end,3); % Sampling frequency downsamplingFctr = 1; Ts = mean(diff(t_raw))*downsamplingFctr; % seconds; fs = 1/Ts; [x,t] = resample(x_raw,t_raw,fs); y = resample(y_raw,t_raw,fs); % Normalize output y = y/sum(y)*sum(x); figure; plot(t,x); hold on; plot(t,y); xlabel('time [s]'); X = fft(x); Y = fft(y); H = Y./X; h = ifft(H); % Define discrete frequency axis f_axis = [-length(X)/2:length(X)/2-1]/length(X)*fs; figure; plot(f_axis,fftshift(abs(X))); hold on; plot(f_axis,fftshift(abs(Y))); plot(f_axis,fftshift(abs(H))); xlabel('frequency [Hz]'); ylabel('Magnitude'); legend('Input','Output','Estimated transfer function') figure; plot(real(h)); hold on; xlabel('n') ylabel('Impulse response h[n]');

enter image description here

In the case of the first signal set, things are trickier because the duration of the recoded signal is much longer. If you use the same approach as above, the bandwidth of the input signal becomes much smaller than the sampling frequency, which results in a terribly noisy estimate. Therefore I played around with downsampling factor downsamplingFctr. If you set values > 100 you start seeing a sensible impulse response, but you still have a significant amount of noise.

Note that downsampling implies low-pass filtering of the signals. Therefore, you can't just increment the downsampling factor indefinitely, because you would end up with significant distortion of the input signal and the result would be incorrect, even if it "looks smoother".

All in all, system identification through simple FFT division has its drawbacks and is in general deprecated as far as I know (there are specific and motivated exceptions). You will probably get better results with other techniques (Wiener filtering might be a good keyword).

Improve this answer
$\endgroup$
6
  • $\begingroup$ Ha, good catch on the division. $\endgroup$ Commented Nov 3, 2022 at 21:07
  • $\begingroup$ Hi Vito, all the values, also the original data are descrete values. So if I speak about integration I of course mean discrete Integration/summation. But to calculate the area under the RTD-curve you have to consider the values of the time steps dt or delta t also. $\endgroup$ Commented Nov 3, 2022 at 21:11
  • $\begingroup$ I'm not sure if I need an element-wise divison or if the devonvolution is performed element-wise. The resulting RTD looks better without the element-wise division. $\endgroup$ Commented Nov 3, 2022 at 21:13
  • $\begingroup$ If I performe the divison COUT/CIN this results in a n by n matrix for RTD, where n equals the length of the vectors COUT or CIN. In the first column of the RTD-matrix I find the desired values. The other columns are all zero. Therefore rtd=rtd(:,1). $\endgroup$ Commented Nov 3, 2022 at 21:18
  • 1
    $\begingroup$ Regarding the division, this is what I was fearing. You are probably dividing a column by a row, not sure how MATLAB interprets that. Maybe you end up "by chance" with a correct result by retaining the first column. But you definitely need an element-wise division. $\endgroup$ Commented Nov 3, 2022 at 21:30
1
$\begingroup$

I postet the question as a guest and proceeded with creating an account with the same email address after that, but unfortunately it didn't connect this account with the guest account. So I'm not able to comment or edit for the moment.

This is a response to @Vitos comment on the elementwise division under answer 1.

If you divide a vector by another vector in matlab (A/B, not element-wise), it solves a system of linear equations for $x$ to satisfy $A= xB$. (see documentation for the matlab function mrdivide).

When using this operation somehow my results turn out good and make sense. But when I use the element-wise division A./B they don't make sense.

And this is a response to @Vitos comment on the summation

You are right the sum of the deconvolved impulse response $RTD$ (sum(rtd)) is always 1, independent on the sample frequency $f_s$. And if you integrate it with resp. to the time steps determined by $f_s$ it takes up the value of $1/f_s$.

So the deconvolution in the frequency domain COUT/CIN obviously does not account for the time step information. Yet this information is relevant.

If I conduct a conventional deconvolution by solving a system of linear equations I'd also consider the time steps and the resulting rtd would have the unit [1/s].

(For example the first equation of an nummerical convolution would be $$ C_{out}[1]=C_{in}[1] RTD[1] \Delta t $$

and the solution for $RTD[1]$ would be

$$ RTD[1]=\frac{C_{out}[1]}{C_{in}[1]\Delta t} $$

$[RTD] = 1/s$ )

So is there a 'mathematical correct' way to include the time step size information in the deconvolution via COUT/CIN?

Here are two examples of the original signals.

Signals1 with a long period https://ufile.io/139k2r58

Signals2 with a short period https://ufile.io/dctxdkzh

For each file the first column is the time vector time_orig, the second column is the input Cin_orig, and the third column is the output Cout_orig.

Please note that due to the process of measurement the area under Cin_orig isn't equal to Cout_orig. So you first have to normalize one signal with the other.

EDIT

OK, I tried it with cpsd pwelch and tfestimate but the results were just zero. I also tried applying a Wiener deconvolution by calculating the transfer function $G(f)$ with $$ \frac{1}{G(f)}=\frac{COUT(f)}{CIN(f)}$$ and then backcalculating the impulse response $H(f)=RTD(f)$ by solving the problem $$ 0=H(f)\biggl(1+\frac{1}{|H(f)|^2SNR}\biggr)-\frac{1}{G(f)}$$ but that gave me worse results than simply using $RTD(f)=\frac{COUT(f)}{CIN(f)}$. So I'll stick to the later approach. I think I can live with the noisy results, if necessary I can smoothen them artificially.

To summarize your answers for my original question: By transforming the signals into the frequency domain and performing a deconvolution, the time step information becomes irrelevant. So the reverse transformed impulse response will give us $h(n)$ with $n$ as an integer control variable. Due to the algorythm the values of $h(n)$ will be scaled that $\sum h(n)=1$ is always valid. For a sample frequency of 1/s this wont make any difference but for a sample frequency of 4/s we have 4 times the ammount of the original values, the original values therefore are scaled down by 4. So to get the time-continous impulse response function $h(t)$ I have to rescale it with the sample frequency: $h(t)=h(n)f_s$. And thereby its physical unit $[1/s]$ also will be restored.

Maybe you guys can confirm this information.

Improve this answer
$\endgroup$
17
  • $\begingroup$ Now I'm confused. sum(rtd) should not always be 1 regardless of the sampling rate, since you're discarding/adding samples by downsampling/interpolating. Unless you're normalizing by length(rtd). As far as the element-wise division, Vito is correct: that's how you should be doing that operation. $\endgroup$ Commented Nov 4, 2022 at 11:52
  • $\begingroup$ I don't think he is changing the sample rate. $\endgroup$ Commented Nov 4, 2022 at 11:54
  • $\begingroup$ "When using this operation somehow my results turn out good and make sense": if you're not sure why you're getting the result you want, there's a problem. Maybe you're getting something that seem to make sense visually, but that is coming from faulty operations, and tricking you into thinking you have the correct result. Try the method I mentioned in an earlier comment on your OP and see what you get for comparison purposes. Again, matrix division in the frequency domain does not make sense for what you're trying to achieve. $\endgroup$ Commented Nov 4, 2022 at 11:56
  • 1
    $\begingroup$ That's difficult to answer without a knowledge of the matter. If your sampling frequency is too low, you might be violating Nyquist's theorem. If it is too high, you might be picking too much noise. $\endgroup$ Commented Nov 4, 2022 at 13:10
  • 1
    $\begingroup$ @Jdip I attached some .csv files of the original data on my answer $\endgroup$ Commented Nov 4, 2022 at 14:23

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.