On Wednesday, May 8, 2013 10:53:09 PM UTC+12, Francesco Perrone wrote:
> The final goal I am trying to achieve is the generation of a ten minutes time series: to achieve this I have to perform an FFT operation, and it's the point I have been stumbling upon.
>
>
>
> Generally the aimed time series will be assigned as the sum of two terms: a steady component U(t) and a fluctuating component u'(t). That is
>
>
>
> u(t) = U(t) + u'(t);
>
>
>
> So generally, my code follows this procedure:
>
>
>
> 1) Given data
>
>
>
> time = 600 [s];
>
>
>
> Nfft = 4096;
>
>
>
> L = 340.2 [m];
>
>
>
> U = 10 [m/s];
>
>
>
> df = 1/600 = 0.00167 Hz;
>
>
>
> fn = Nfft/(2*time) = 3.4133 Hz;
>
>
>
> This means that my frequency array should be laid out as follows:
>
>
>
> f = (fn+df):df:fn;
>
>
>
> But, instead of using the whole f array, I am only making use of the positive half:
>
>
>
> fpos = df:fn = 0.00167:3.4133 Hz;
>
>
>
> 2) Spectrum Definition
>
>
>
> I define a certain spectrum shape, applying the following relationship
>
>
>
> Su = (6*L*U)./((1 + 6.*fpos.*(L/U)).^(5/3));
>
>
>
> 3) Random phase generation
>
>
>
> I, then, have to generate a set of complex samples with a determined distribution: in my case, the random phase will approach a standard Gaussian distribution (mu = 0, sigma = 1).
>
>
>
> In MATLAB I call
>
>
>
> nn = complex(normrnd(0,1,Nfft/2),normrnd(0,1,Nfft/2));
>
>
>
> 4) Apply random phase
>
>
>
> To apply the random phase, I just do this
>
>
>
> Hu = Su*nn;
>
>
>
> At this point start my pains!
>
>
>
> So far, I only generated Nfft/2 = 2048 complex samples accounting for the fpos content. Therefore, the content accounting for the negative half of f is still missing. To overcome this issue, I was thinking to merge the real and imaginary part of Hu, in order to get a signal Huu with Nfft = 4096 samples and with all real values.
>
>
>
> But, by using this merging process, the 0th frequency order would not be represented, since the imaginary part of Hu is defined for fpos.
>
>
>
> Thus, how to account for the 0th order by keeping a procedure as the one I have been proposing so far?
The process is called "phase randomisation" (or randomization) and there have been several posts on it in in the newsgroup. It is widely used in Monte Carlo modelling, especially of turbulence, where it is used to determine the extent to which an actual signal contains structures. Phase randomisation destroys all structures, but retains the 2nd order properties of a signal.
It makes no sense to use a normal distribution for phase  it is a circular quantity: use uniform instead.
phi=rand(Nfft/2,1)*2*pi;
nn=exp(i*phi);
Hu=Su.*nn;
To form the Fourier transform, whence the time series, do this:
Y=[0;Hu;flipud(conj(Hu))];
y=ifft(Y);
Note that this gives Nfft+1 data points.
Clearly, if you want exactly Nfft points, you must adjust your algorithm accordingly.
Note also that if Su is a true PSD, then you must transform into the amplitude spectrum by sqrt(Su*df)
