Apply Phase Shifts Using FFT for Gaussian Noise Signal

[marksayles]

What I want to do is apply a frequency-independent phase shift to a Gaussian Noise signal. Obviously I can just invert it to get a pi shift. Also I can take the Hilbert transform to get a pi/2 shift and I can invert that to give myself in effect a 3/2pi shift. BUT... what I want to do is to be able to apply any phase shift, 10, 20, 30 degrees etc.

The way I have been thinking about this is that I can take the FFT of the noise and get the magnitude and phase information, and then I could manipulate the phase but leave the magnitude the same and take the IFFT to get my phase-shifted signal. Here is an example bit of MATLAB code doing this on a single sinusoid. As it is here the code doesn't work. But if you change the phase shift p from pi/2 to pi then it does? Is there something I'm missing. Perhaps I need to apply the phase shift differently to the positive and negative parts of the FFT? Or is it something to do with the fact that theta is distributed on [-pi pi] and when I add "p" I change this and make it ambiguous or something?

sr=10000;
dt=1/sr;
len=0.5;
t=0:dt:len;
f=500;
%generate signal
sig=sin(2*pi*f*t);
%Define a phase shift in rads
p=pi/2;
%Get the FFT of the signal
z=fft(sig);
%Get the magnitude
mag=abs(z);
%Get the phase
theta=angle(z);
%Add the phase shift onto the phase
theta=theta+p;
%Set the new real and imag parts of the spectrum
R=mag.*cos(theta);
I=mag.*sin(theta);
%Make the new spectrum
spec=complex(R,I);
%Get the new signal
newsig=real(ifft(spec));
%plot the signals
figure;plot(t,sig);hold on; plot(t,newsig,'g');

 

[dimensionless]

Why doesn't your code work?

May I ask where the following lines come from?

R=mag.*cos(theta);
I=mag.*sin(theta);

 

[marksayles]

The code works if the phase shift (p) is set to pi, but not for other values. Which makes me think I've made some goof with the maths and/or my fundamental understanding of this topic, rather than a simple MATLAB slip.

The two lines come from the fact that I want to keep the magnitude the same but just change the phase. So the magnitude of the complex number z = a +bi is sqrt(a^2+b^2). I calculate the angle of each component in theta and then add the phase shift in rads. Viewing the complex number spec on polar coords we have an x-axis which is cos(theta)*mag and a y-axis which is sin(theta)*mag.

I think this is right... but would be very happy to hear why I'm wrong if that is the case.

 

[dimensionless]

I haven't run your code, but it looks okay. You might try doing an N-point fft, where N is

N = 2^nextpow2(length(sig));

and then

z=fft(sig, N);

and then later on

newsig=real(ifft(spec, N));

I don't know what N defaults too or what effect the default might have. The fft (and mag), though, should produce a single spike at the frequency of sig. Does newsig have the same frequency as sig?

 

[SammySmith]

frequency independent phase shift

hey marksayles,

I have a similar problem as you..Have you figured out how to get 'any phase shift' using FFT and then IFFT.
I will be waiting for your reply.

Sammy

 

[mharaghi]

the problem arises because the signal loses it symmetry, and the ifft gives a complex answer. it is obvious that a shifted sine doesn't have a complex component...

You can test it, by manually adding a delay to the inpt sine wave and see it's symmetric fourier.


but by adding the theta, it will lose it's symmetry. I don't know how to solve it, but that's the problem.

 

[mharaghi]

this code is the way it should be implemented, (a freq component was missed in your code).

This is still not correct, and I am missing something, some coefficient ...

clc
sr=10000;
dt=1/sr;
len=0.05;
t=0:dt:len;
f=500;
%generate signal
sig=sin(2*pi*f*t);
%Define a phase shift in rads
p=-pi;
%Get the FFT of the signal
z=fft(sig);
%Get the magnitude
for k=1:f+1
spec(k)=z(k)*exp(-i*2*pi*k*p);
end

%Get the new signal
newsig=(ifft(spec));
%plot the signals
figure;plot(t,sig);hold on; plot(t,newsig,'g');

 

[ldvbin]

See code below, a couple of notes:
Phase shift with a FFT can be looked at from a sample delay/FFT perspective. Using the second and third equations from the following table (compare to the Z-transform sample delay):
http://en.wikipedia.org/wiki/Discrete_Fourier_transform#Some_discrete_Fourier_transform_pairs
We can delay the original sine wave by a give number of samples. All that is needed then is to figure out how many samples comprise the desired phase shift. Another quick note is that the number of samples shifted needs to be an integer or the amplitude of the sine wave is affected due to the rectangular windowing effects inherent in the DFT (FFT).

clc
clear
sr=10000;
dt=1/sr;
len=0.01;
t=0:dt:(len-dt);
f=500;
N = length(t);

%generate signal
sig=sin(2*pi*f*t);

%Define a phase shift in rads
p=-pi/4;
num_samp = round((sr/f)*(p/(2*pi)));

%Get the FFT of the signal
z=fft(sig);

%Delay each fft component
for k=1:length(z)
w = 2*pi/N*(k-1);
spec(k)=z(k)*exp(-j*w*num_samp);
end

%Get the new signal
newsig=(ifft(spec));

%plot the signals
figure;plot(t,sig);hold on; plot(t,newsig,'g');

 

[mea03wjb]

Please can you define what the varaibles are at the top of this script:

What are 'sr' and 'f'?

 

[astrogirl123]

I have a questions regarding the same problem :

I have a spectrum with flux on the Y-axis and wavelength on the X-axis. What I want to do is take the Fourier transform of this spectrum. Then add a random phase between 0 and 2pi to the phase only. Then take the inverse Fourier Transform of this. The piece of code I wrote in IDL for a simple example is below. I was hoping to scramble up the spectrum by doing this. But I am not sure if I have achieved that goal. I have two questions

1. The final spectrum ( for some random phases) are vertically offset from the original spectrum. If the DC component is zero in the original spectrum, then this problem doesn't exist! Why is that so?

2. The amplitude of the final signal varies a lot from the original spectrum. I wasnt expecting that either. why is that?

3. All I want is that the components of the original spectrum be jumbled up in the X-direction (ie the wavelength direction). How can it be done?

MY CODE :

n = 256
x = FINDGEN(n)
y = COS(x*!PI/6)*EXP(-((x - n/2)/30)^2/2)+1.00

tek_color

yfft = fft(y)

magnitude = abs(yfft)
angle = ph(yfft)

for count=0, 300 do begin

rand = randomu(seed, 1)*2*!pi
randph = replicate(rand, 256)

fft_signal = magnitude*exp(complex(0,1)*(angle+randph))
ifft_signal = (fft(fft_signal, /inverse))

wait, 0.2

plot, x, y
oplot, x, (ifft_signal), color=2

endfor

end

 

[AlephZero]

1. The final spectrum ( for some random phases) are vertically offset from the original spectrum. If the DC component is zero in the original spectrum, then this problem doesn't exist! Why is that so?
end

Read the documentation on exactly how IDL stores the fft data.

I don't have IDL so I don't know the answer to that question, but the general issue is that if you do an FFT on N data points where N is even, there are N+1 frequency components in the FFT not N. But the lowest (0) and highest (the Nyquist frequency) only have an "in phase" (or real) component and not a quadrature (or imaginary) component.

So, many FFT routines store the real part of the N+1'th term where the imaginary part of the 0'th term should logically go. That way the FFT data fits into an array of length N, the same length as the input data.

Once you know how IDL handles that, you then need to think about what (if anything) it means to "do a random phase shift" on these two freqency components.

I expect that will explain where your "random level shifts" are coming from, if you are phase-shifting a large zero component. It might also explain your question 2.

 

[astrogirl123]

Thanks for the reply!. What I am looking for is the solution posted by "ldvbin". That code works in MATLAB but for some reason the exact algorithm doesn't work in IDL. In either case, can some one help me with how to relay that information to a data where I have flux data values against wavelength data values. The Y-axis is not defined by a simple functional form of the X-axis. In particular , how to apply the sample delay mentioned by "ldvbin"?