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.
> 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.
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)
> 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
> 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?
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.
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 need to multiply bin 2 by your complex shift variable and bin 255 by the complex conjugate of your shift variable ( just invert the imaginary part). Then when you take the ifft the result should be real only (within numerical precision limits).
>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).
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)
> 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).
Bob, do you see what this MATLAB convention has done to us? now *we* are using it, but because we know it's wrong, we also have to add a note, a caveat or qualification so that people know for sure what we are talking about. either that or risk being misunderstood by 1.
>>> because we know it's wrong
> Why is it wrong?
> There isn't a single known madhouse in history where the numbering of
> padded cells started with 0 ...
I agree the numbering system is bad. I guess it was decided a long time ago and then no one had the courage to break all the existing scripts out there.
> I agree the numbering system is bad. I guess it was decided a long time ago and then no one had the courage to break all the existing scripts out there.
long ago (sometime in the 90s, when MATLAB was less than a decade old) i showed Cleve exactly how to do it and be perfectly backward compatible. we posted it here on comp.dsp and on comp.soft-sys.matlab and there was a lot of discussion. i think Eric Jacobson might remember, but i cannot remember who else was around back then who are still around now. maybe Jerry Avins.
anyway, inside each MATLAB variable (which is an N-dimensional array), there exists now a vector of length N that defines the length of each dimension in the array. the total number of elements in the MATLAB variable is the product of all of the elements of this vector. you can get the values of this vector with the size() function and change the values of this vector with the reshape() function (as long as the product is the same) without touching the elements themselves.
what i proposed was another vector of the same size that defines the index of the very first element of each dimension in the array. the default values when a MATLAB variable is created are always equal to 1. THAT MAKES IT PERFECTLY BACKWARD COMPATIBLE. so no existing scripts are broken. then, to make use of this would be two new functions that i would call "origin()" and "reorigin()" or "base()" and "rebase()". the first would be like size() and tell you what the origins are for each dimension and the latter would be like reshape() and allows you to modify the origins away from their value which has default of 1.
of course some existing math functions would be rewritten to make it useful to any arrays of modified origin. to add or subtract arrays, not only must the size() be the same, also the origin() would have to match or a run-time error occurs. to multiply arrays there would have to be an index match where the column index of one multiplicand array must match the row index of the other.
doesn't require any courage at all. just a little bit of corporate will at the highest level (like Cleve) and a measure of good sense.
>>> because we know it's wrong
>Why is it wrong?
>There isn't a single known madhouse in history where the numbering of
>padded cells started with 0 ...
Numbering the cells from zero is unusual, but numbering the floors from
zero is quite common in many countries. Maybe the populations in those
place are more numerate.
On Mon, 03 Sep 2012 04:36:18 -0500, Lightbulb85 wrote:
> 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:(N-1))/N);
> and found the same results as my earlier code... ie, a scaled version of
> the input, not a phase shifted one.
What I suggested and what Bob Adams told you to do should be exactly the same thing in the case where you only have two non-zero bins -- did you apply the shift I suggested to _all_ the bins?
Also, if you only have exactly one frequency present then the easiest way to phase shift it is to just delay it by a known amount (or if integer delays isn't good enough, some filter of the form (a * z + b) / z^q, with a, b and q chosen for unity gain and the desired shift at the frequency of interest).
-- Tim Wescott
Control system and signal processing consulting
www.wescottdesign.com
