0
$\begingroup$

I have a function with an equation:

$$C = 1.6925\left( e^{-0.136t}-e^{-1.192t}\right) $$

Where $C$ is real and $t$ represents time in hours. Beneath is the representation of my function.

Represents a pharmaceutical drug concentration increase and decay in the body from absorption point at t=0

I am trying to find the best sampling step (based on Nyquist-Shannon theorem) to sample and reconstruct this function using a Fourier Transform.

After using a Fast Fourier Transform algorithm in MATLAB (fft), by setting:

x= fft(c);
xmag=abs(x); (amplitude estimation)

I find the below spectrum for the function for a $t$ between $0$ and $100$ sampled at $1000 \text{Hz}$ , with the highest amplitude at:

bin(1) = 1.1025e+04

and the second highest at bin(100001) with an amplitude of 9.9946e+03. the bin number 100001 represents also the length of xmag.

So my question is how would you interpret the spectrum? Specifically I would like to know how to determine the highest frequency at which I can sample this function and reconstruct it without aliasing? enter image description here

Improve this question
$\endgroup$
  • $\begingroup$ It's not about spectrogram but spectrum. $\endgroup$ – jojek♦ Jun 11 '15 at 19:17
  • $\begingroup$ Yes, you are right thanks also for editing my formulas i did not know how to do that. $\endgroup$ – Wazaa Jun 11 '15 at 19:44
2
$\begingroup$

First, since $t>0$ in your case, you can write your function as

\begin{equation} C(t)=1.6925\left(\exp^{-0.136t}- \exp^{-1.192t} \right) u(t) \end{equation} where $u(t)$ is a unit-step function. Then, denote the FT of $C(t)$ as $C(F)$, which is given by \begin{equation} C(f)=1.6925\left(\frac{1}{0.136+2\pi jf}- \frac{1}{1.192+ 2\pi jf} \right) \end{equation}

The spectrum $C(f)$ is plotted below. We can see that $C(f)$ is almost zero after 5 Hz. You can choose the sampling frequency according to the Nyquist criterion as $f_{\max} \geq 10$ Hz. You have taken a very high sampling rate (1000 Hz) that is 100 times more than the required rate. When I took 10 Hz to sample the signal above and took the FFT I got the results as follows (in order time domain $C(t)$, frequency domain signal $C(f)$, frequency domain domain signal obtained via FFT, that is, FFT(C(t)). Note the inherent periodicity of the FFT spectrum)

Fig.1 Fig. 2 Fig. 3

And here is the MATLAB/Octave code for your reference..

clear all; clc;close all;

%% Time-domain signal
Ts=0.1; % Sampling time (obtained from the frequency domain plot)
t=0:Ts:60; % in sec
C=1.6925* (exp(-0.136*t)-exp(-1.192*t));
figure,plot(t,C);
title('Time domain signal');
xlabel('Time in sec')

%Frequency domain signal
f=0:.1:10;
F1=1.6925*((1./(0.136+(j*2*pi*f)))-(1./(1.192+(j*2*pi*f))));
figure,plot(f,abs(F1));
title('Frequency domain signal');
xlabel('Frequency in Hz')
% Via FFT
F2=fft(C);
L=length(F2);
figure,plot((0:L-1)*(1/(L*Ts)),abs(F2));
title('Frequency domain signal via FFT of C(t)');
xlabel('Frequency in Hz')
Improve this answer
$\endgroup$
  • $\begingroup$ Could you please clarify your explanation, my time is in hours but yours is in seconds but we have almost the same graph and also how are you able to determine the frequency when FFT in Matlab only gives you back bin numbers and amplitudes while in your frequency representation your abscissa is in hertz? $\endgroup$ – Wazaa Jun 12 '15 at 16:39
  • $\begingroup$ The frequency bins are converted into frequencies with the formula :frequency = (sampling frequency)*(frequency bin number)/(total number of samples) $\endgroup$ – Seetha Rama Raju Sanapala Jun 12 '15 at 23:49
  • $\begingroup$ Yes Wazaa, Seetha Rama Raju is correct..It is just a scaling aspect very common in signal processing...Each bin in the FFT domain corresponds to normalized digital frequency (normalized to sampling frequency). To find the un-normalized continuous frequency we do the transformation mentioned in the code. $\endgroup$ – Oliver Jun 13 '15 at 0:53
  • $\begingroup$ Regarding the time in hrs question, it is just a matter of scaling by 3600 sec = 1 hr..Let me know if you cannot figure it out..Hint: Start with the continuous-time Fourier transform integral and note down the units of time and frequency $\endgroup$ – Oliver Jun 13 '15 at 0:55
  • $\begingroup$ Thank you once again guys for your explanations.....Questions: 1. I don't find in the litterature (MIT Oppenheim's book) a reference to the fact that the highest frequency is the one at which C(f) is zero, could you point me to one? 2. I see in the C(f) graph that C(f) is nullified even earlier (between 1 and 2), yet both of you are referring to 5 Hz as the zero point so I am guessing here that i must have missed something, what is ?. $\endgroup$ – Wazaa Jun 29 '15 at 15:55
-1
$\begingroup$

How did you arrive at sampling frequency of 1000 Hz? Ideally speaking, you should find the Fourier transform of the given function analytically. This can be done using the FT tables for exponentials. Then you should find the highest frequency beyond which the FT is almost zero i.e less than the smallest number represented by finite precision you are using. Twice or more of this frequency is your sampling rate. How many samples should you take? The signal is almost zero beyond 5 x 10$^4$ as can be seen from the figure. So you take samples upto this time point. Now you take the FFT of these samples and it would be equal to the FT you analytically calculated i.e it would be sampled version of the FT. If you used sampling rate less than this procedure explained, the resulting FFT would be aliased version.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.