I did not understand everything, because I do not know R. Can you explain it without pseudo-code?

You can send a mail [email protected]

Here is how I would do it in Matlab:

Generate a signal that consists of two sinusoids of frequencies 15 Hz and 40 Hz. The signal is sampled at 100 Hz for one second.

fs = 100; t = 0:1/fs:1-1/fs; x = sin(2*pi*15*t) + sin(2*pi*40*t);

Compute the discrete Fourier transform of the signal. Plot the magnitude of the transform as a function of frequency.

y = fft(x); ly = length(y); f = (-ly/2:ly/2-1)/ly*fs;

Compute the phase:

phs = angle(fftshift(y));

Compute the phase of the peak amplitude fft:

[~,maxindex]=max(abs(y));

maxphase=phs(maxindex);

Generate random phases:

phasemin = 0; phasemax = 2*pi; fftlength=length(y);

newphase = (b-a).*rand(fftlength,1) + a;

Inject the previous phase of the maximum peak in fft:

newphase(maxindex)=maxphase;

Note that this is for a double sided spectrum. The code needs some adapting, seeing as it is done for a double sided spectrum, but the approach is similar. Hope this helps.

Below I describe my solution that I sent already to Iman, but that can be useful for other people

The solution is not very complicated.

1 - Assume you have computed the Fourier transform defined by the vector X(k), k=0,1 ..., N-1.

2 - For simplicity assume that N is even and that the frequencies you want to retain are 3 at the subscripts n_p =[n1 n2 n3]. Let

M=[0, n1, n2, n3, N/2]

3 - Let P be the vector with the phases P(k) for k=0, 1, ..., N/2.

4 - Create a new vector, Q, with the phases

[P(0)=0, P(n1), P(n2), P(n3), P(N/2)=0 or pi ].

5 -  Remove the phases in Q, from P, starting from right to left:

the first to remove is the phase at N/2, then n3 ... To remove a given phase at "j" do P=[P(0), ..., P(j-1), P(j+1), ...]

6 - Produce the randomization of the remaining P.

7 - Reinsert the previously removed phases starting at 0 and continuing to right. To insert  a given phase at "j" do P=[P(0), ..., P(j-1), Q(n_i), P(j), ...]

8 - Reconstruct the total phase vector having in mind the hermiticity:

P(N/2+k) = conj(P(N/2-k), k=1, 2, ..., N/2-1

Do not trust much on the "spectrum" for computing the frequencies.

Dear @Costin Vasile and Dear @Manuel Duarte Ortigueira,

Thanks a lot for your nice answers. I was actually increasing my personal knowledge about Fourier Transform and complex numbers as the fundamentals to apply your approaches.

So I practically applied some of your tricks and also some theories as follow. Maybe it can be hopefully helpful to others facing the same problem, or if they can improve or revise it;

If 'Z' is a complex number which is generated by applying Fast Fourier Transform: fft(n) = Z

Z = a + bi in which 'a' is the Real part and 'b' is the Imaginary part of the complex number 'Z'.

If we call 'r' as the magnitude or Modulus of 'Z', then it can be calculated as the absolute value of 'Z' or;

r = |z| or r = √a^2+b^2

And if we call '∂' as the argument of 'Z' or the angle between 'a' value on the X-axis and 'b' on the Y axis which produces the vector of 'Z', then we can use the following equations ;

∂= arctan(b/a)

a = r cos ∂

b = r sin ∂

Z = r (cos ∂ + i sin ∂)

or in exponential form;

Z = r e^i ∂

Then if we are going to randomize phases or in the other word the angle between 'a' and 'b' for each FFT component the ∂ will be replaced by another random value between 0-2π, and we call it ∂', so;

Z' = r cos ∂' + r sin ∂' i

for which the 'r' can be calculated using { r = √a^2+b^2 } equation.

So I used a for() loop to randomize the phases for the length of all FFT components by generating a random number between 0˚-360˚ or 0-2π for each component. I wrote the code in R as follow;

{R}

out