Phase shift in the frequency domain - Problems!!

Hi,
I apologise for the rather basic question.

I am attempting to use MATLAB to carry out a phase shift in the frequency
domain by the following process:

I have a time domain signal, x, with the following parameters:
fs =256;             % Sample Frequency
N = 256;             % Size of FFT
l = 256;             % Length of data

t=[0:1/fs:l/fs-1/fs]'; % Time vector is 256 points long

x=cos(2*pi*2*t+pi/4); % 2 Hz cosine wave with 45 degree phase shift

I compute the FFT of the input signal 
X=fft(x,N)/N;
and get 0.3536+0.3536i in the 2 Hz bin as expected.

I understand that a phase shift can be achieved in the frequency domain by
multiplying my complex X values by cos(P)+jsin(P) where P is my target
phase shift.

Therefore I generate:
shift=complex(cos(pi/4),sin(pi/4));
which gives me
shift = 0.7071 + 0.7071j;

I then multiply X by shift as follows:
X_SHIFTED=X*shift

when I check the following I find
angle(X(3)) = 0.7854          % The angle of the 2 Hz bin of the input
and
angle(X_SHIFTED(3))=1.5708    % The angle of the 2 Hz bin of the shifted
data.

This indicates that I have successfully achieved my target pi/4 phase
shift.  However it is here that I run into problems.

I want to plot the ifft of X_SHIFTED along with my original input to
compare and visualise the phase shift.

When I compute
output=ifft(X_SHIFTED)*N;
and plot, I expect output to be the same magnitude as x but shifted by 45
degrees.  However I find that it is actually a scaled version of x with the
same phase.

I am not sure what I am doing wrong because everything seems to work OK up
until I compute and plot the IFFT.  I am sure it is some glaringly obvious
thing I  have missed.

Any help would be much appreciated.

Am 31.08.12 15:27, schrieb Lightbulb85:
> Hi,
> I apologise for the rather basic question.
>
> I am attempting to use MATLAB to carry out a phase shift in the frequency
> domain by the following process:
>
>
> x=cos(2*pi*2*t+pi/4); % 2 Hz cosine wave with 45 degree phase
> X=fft(x,N)/N;
> shift=complex(cos(pi/4),sin(pi/4));
> shift = 0.7071 + 0.7071j;
> X_SHIFTED=X*shift
> output=ifft(X_SHIFTED)*N;
> and plot, I expect output to be the same magnitude as x but shifted by 45
> degrees.  However I find that it is actually a scaled version of x with the
> same phase.

output is not anymore real, it is a complex vector. Are you plotting 
real and imaginary parts? Magnitude? If you want pure real output, you 
must treat the negative frequencies differently from the positive one, 
such that they are remain complex conjugates. Try it by comparing the 
FFT of your expected output with the product X_SHIFTED.

	Christian

If you use a complex input: 

x=complex(cos(2*pi*2*t+pi/4),sin(2*pi*2*t+pi/4));

you'll get the results you're looking for.

On Friday, August 31, 2012 11:38:07 AM UTC-4, David Drumm wrote:
> If you use a complex input: 
> 
> 
> 
> x=complex(cos(2*pi*2*t+pi/4),sin(2*pi*2*t+pi/4));
> 
> 
> 
> you'll get the results you're looking for.

Lightbulb:
The way you are computing your shift is wrong. The shift is frequency dependent. As far as I can remember, if you want x(t-t0), you need to multiply by exp(-j 2pi f t0) 

- Fernando

On Friday, August 31, 2012 12:00:29 PM UTC-4, f.qu...@gmail.com wrote:
> On Friday, August 31, 2012 11:38:07 AM UTC-4, David Drumm wrote:
> 
> > If you use a complex input: 
> 
> > 
> 
> > 
> 
> > 
> 
> > x=complex(cos(2*pi*2*t+pi/4),sin(2*pi*2*t+pi/4));
> 
> > 
> 
> > 
> 
> > 
> 
> > you'll get the results you're looking for.
> 
> 
> 
> Lightbulb:
> 
> The way you are computing your shift is wrong. The shift is frequency dependent. As far as I can remember, if you want x(t-t0), you need to multiply by exp(-j 2pi f t0) 
> 
> 
> 
> - Fernando

I forgot to add: if you want x(t-t0), you need to multiply by exp(-j 2pi f t0) in the frequency domain

On Fri, 31 Aug 2012 08:27:33 -0500, Lightbulb85 wrote:

> Hi,
> I apologise for the rather basic question.
> 
> I am attempting to use MATLAB to carry out a phase shift in the
> frequency domain by the following process:
> 
> I have a time domain signal, x, with the following parameters: fs =256; 
>            % Sample Frequency N = 256;             % Size of FFT
> l = 256;             % Length of data
> 
> t=[0:1/fs:l/fs-1/fs]'; % Time vector is 256 points long
> 
> x=cos(2*pi*2*t+pi/4); % 2 Hz cosine wave with 45 degree phase shift
> 
> I compute the FFT of the input signal X=fft(x,N)/N;
> and get 0.3536+0.3536i in the 2 Hz bin as expected.
> 
> I understand that a phase shift can be achieved in the frequency domain
> by multiplying my complex X values by cos(P)+jsin(P) where P is my
> target phase shift.

Fernando was right: you need to multiply point-by-point with a frequency-
dependent phase shift.  For an N-point FFT and a time shift in samples, 
you need:

x_shifted = x .* exp(2 * pi * t_shift * (0:(N-1)) / N);

Note that this works "perfectly" for integer shifts, and less so for non-
integer shifts (because of the phase discontinuity at high frequencies).  
If you're windowing the original data anyway, though, I'm pretty sure 
that a part-sample shift will work out mostly OK.

Also note that I'm not sure of the sign convention on the t_shift -- 
t_shift positive may shift into the future or past the way I've expressed 
it above.  Check.

-- 
My liberal friends think I'm a conservative kook.
My conservative friends think I'm a liberal kook.
Why am I not happy that they have found common ground?

Tim Wescott, Communications, Control, Circuits & Software
http://www.wescottdesign.com

Hi,
Thanks everyone for your comments, they are very helpful.

I have a few comments of my own:

Fernando and Tim:
Thanks for your comments.  I understand that the phase shift is frequency
dependent.  In my example, there is only one frequency with any energy
(2Hz) so surely it is not necessary to shift every frequency by a different
amount, just as long as I shift the 2 Hz content by my target amount?  Or
am I missing something?  Do I need to shift the negative frequencies by a
different amount?

I implemented the code Tim suggested:
X_SHIFTED=X.*exp(2*pi*t_shift*(0:(N1))/N);

and found the same results as my earlier code... ie, a scaled version of
the input, not a phase  shifted one.

It is also noteworthy that when I checked the angle of the signals, it
seemed to work using both methods.  But when I calculate the ifft and plot,
the phase of the input and output appear to be the same.

Christian:
I think you are right that it is because the output is not real any more. 
When I plot, I get the MATLAB warning that the imag parts are being
ignored.  I am plotting only the real part.  I think the phase shift
worked, it is just how I am attempting to visualise it.  Do you have any
suggestions how I might visualise it correctly?

David:
You are also correct, when I used a complex input I got the results I
wanted.

However, this leads me to another question... how to achieve this in real
life...?

The reason behind my question is that I am trying to further my
understanding of the FFT and in particular, the twiddle factors. My example
above is simply to help me illustrate and understand the shifting
algorithm.

So for the FFT...
I understand that we use bit-reversal to break the data (say an 8-point
time domain signal) up into 8 separate frequency domain signals.

Then the complicated part is reassembling them into one 8-point frequency
domain signal.

So I now have 8 frequency domain points.
F0
F1
F2
F3
F4
F5
F6
F7

I want to combine F1 and F2 so I must stretch the signals.  In the time
domain, this would look like the following:
[t0 0]
[t1 0]

If I wanted to add these in the time domain I would now shift the "odd"
data points and get
[t0 0]
[0  t1] which can now be added to produce my target [t0 t1].

So in the frequency domain I stretch my signals and get
[F0 F0]
[F1 F1]

and now I want to phase shift the "odd" signal to achieve my time domain
time shift.

So I am attempting to better understand phase shifting in the frequency
domain... hence my earlier post.

In real life my input will be completely real.  
[t0 t1 t2 t3 t4 t5 t6 t7]

David suggested I must use a complex input.  So in real life I believe I
need to insert zeros into the imaginary part.

I am struggling to bridge the gap between theory and practice.  So I now
have two vectors:
[t0 t1 t2 t3 t4 t5 t6 t7] and [0 0 0 0 0 0 0 0]
How do I now actually go about beginning the FFT...?  I can bit-reverse
[t0...t7] but what happens to my imag part (ie, the zeros)?
Once I have 8 separate frequency domain signals, I can phase shift using
the previously posted method suggested by Tim and Fernando or by David
(even though I am still struggling to visualise the effect in MATLAB).  But
I’m not sure how to get started.

Thanks for everyone's suggestions so far.

If you examine your fft output you will find that there are 2 bins that are=
 non zero, bin 3 and bin 255 (assuming the bins run from 1 to 256). You nee=
d to multiply bin 2 by your complex shift variable and bin 255 by the compl=
ex conjugate of your shift variable ( just invert the imaginary part). Then=
 when you take the ifft the result should be real only (within numerical pr=
ecision limits).

Bob

>If you examine your fft output you will find that there are 2 bins that
are=
> non zero, bin 3 and bin 255 (assuming the bins run from 1 to 256). You
nee=
>d to multiply bin 2 by your complex shift variable and bin 255 by the
compl=
>ex conjugate of your shift variable ( just invert the imaginary part).
Then=
> when you take the ifft the result should be real only (within numerical
pr=
>ecision limits).
>
>Bob
>
>
>
Hi Bob,
Done and works, thank you so much.

Hi,

some copy&paste delay code 
- waveform is cyclic
- rate_Hz: sampling rate (redundant, use 1)
- delay_s: desired delay in s (or samples if using rate_Hz = 1)

used here:
http://www.dsprelated.com/showcode/288.php

Cheers

Markus

% ****************************************************************
% delay by phase shift
% needed only for demo code
% ****************************************************************
function waveform = cs_delay(waveform, rate_Hz, delay_s)
    rflag = isreal(waveform);
    
    n = numel(waveform);
    cycLen_s = n / rate_Hz;
    nCyc = delay_s / cycLen_s();
    
    f = 0:(n - 1);
    f = f + floor(n / 2);
    f = mod(f, n);
    f = f - floor(n / 2);
    
    phase = -2 * pi * f * nCyc;
    rot = exp(1i*phase);
    
    waveform = ifft(fft(waveform) .* rot);
    if rflag
        waveform = real(waveform);
    end
end