Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
randomize phase of Fourier transform and get real data from inverse Fourier

Subject: randomize phase of Fourier transform and get real data from inverse Fourier

From: Diederick C. Niehorster

Date: 13 Jan, 2011 08:39:50

Message: 1 of 7

Dear all,

I'm puzzling on the following problem:
I want to randomize the phase spectrum of the Fourier transform of an
image, but leave the power spectrum untouched, to get a randomized
image with the same power spectrum as my reference image.

I know that for the inverse transformed to be real, there has to be a
certain symmetry to the Fourier spectrum, as I use a random phase, I
need to recreate this. (Is conjugate symmetry the right word?).

I appear to fail however, with a significant part of my image ending
up in the imaginary part (at least its pretty random, on one run
almost all is imaginary, at other times its pretty close to all real).
So I'm probably recreating something wrongly.

Could you tell me where my thinking is in error so I can correct (and
learn!)?
Thank you very much for your time,
Best,
Dee

% get image from matlab that we all have, see
% http://blogs.mathworks.com/steve/2006/10/17/the-story-behind-the-matlab-default-image/
defimage = pow2(get(0,'DefaultImageCData'),47);
highbit = 51;
lowbit = 47;
numbits = highbit - lowbit + 1;
b = bitshift(a,-lowbit);
b = fix(b);
b = bitand(b,bitcmp(0,numbits));
theimage = b/max(b(:));

%% here's the start of my code:
% do fft
fftimg = fft2(theimage);
amp = abs(fftimg);
ph = angle(fftimg);

% random phase
rph = 2*pi*rand(size(theimage))-pi;
% make conjugate symmetric
%
% AAAAAAAAAAAAAAAAA
% B
% B C
% B
% B F dc E
% B
% B D
% B

phasemat = fftshift(rph);

S = size(rph,1);

A = rph(1,:);
B = rph(2:S,1);
C = rph(2:S/2,2:end);
D = -flipud(C);
E = rph(1+S/2,(2+S/2):S);
F = -fliplr(E);
dc = rph(1+S/2,1+S/2);

srph = ifftshift([A;[B [C;[F dc E];D]]]);

% recreate DFT coeffs, do inverse transform
nfftimg = amp .* exp(1i.*srph);
nimgdata = ifft2(nfftimg);

% test if its all in real
sum(real(nimgdata(:)))
sum(imag(nimgdata(:)))

Subject: randomize phase of Fourier transform and get real data from

From: Diederick C. Niehorster

Date: 16 Jan, 2011 02:17:20

Message: 2 of 7

Dear All,

Sorry for the long post before. As I am no Fourier expert, could you
verify the following to me?

- If i take the ftt of an image and then ifft the result I get a fully
real image again (of course!)
- however, when I first randomize the phase spectrum (taking symmetry
into account, hopefully correctly) and then transforming back, I would
expect my image to still fully lay on the real plane. this appears to
not be true. Is that expected, or is this a coding error on my part?

An answer to that would help me a lot to look further!
Thank you!

Best regards,
Diederick

On Jan 13, 4:39 pm, "Diederick C. Niehorster" <dcni...@gmail.com>
wrote:
> Dear all,
>
> I'm puzzling on the following problem:
> I want to randomize the phase spectrum of the Fourier transform of an
> image, but leave the power spectrum untouched, to get a randomized
> image with the same power spectrum as my reference image.
>
> I know that for the inverse transformed to be real, there has to be a
> certain symmetry to the Fourier spectrum, as I use a random phase, I
> need to recreate this. (Is conjugate symmetry the right word?).
>
> I appear to fail however, with a significant part of my image ending
> up in the imaginary part (at least its pretty random, on one run
> almost all is imaginary, at other times its pretty close to all real).
> So I'm probably recreating something wrongly.
>
> Could you tell me where my thinking is in error so I can correct (and
> learn!)?
> Thank you very much for your time,
> Best,
> Dee
>
> % get image from matlab that we all have, see
> %http://blogs.mathworks.com/steve/2006/10/17/the-story-behind-the-matl...
> defimage = pow2(get(0,'DefaultImageCData'),47);
> highbit = 51;
> lowbit = 47;
> numbits = highbit - lowbit + 1;
> b = bitshift(a,-lowbit);
> b = fix(b);
> b = bitand(b,bitcmp(0,numbits));
> theimage = b/max(b(:));
>
> %% here's the start of my code:
> % do fft
> fftimg  = fft2(theimage);
> amp     = abs(fftimg);
> ph      = angle(fftimg);
>
> % random phase
> rph     = 2*pi*rand(size(theimage))-pi;
> % make conjugate symmetric
> %
> % AAAAAAAAAAAAAAAAA
> % B
> % B       C
> % B
> % B   F   dc   E
> % B
> % B       D
> % B
>
> phasemat = fftshift(rph);
>
> S = size(rph,1);
>
> A = rph(1,:);
> B = rph(2:S,1);
> C = rph(2:S/2,2:end);
> D = -flipud(C);
> E = rph(1+S/2,(2+S/2):S);
> F = -fliplr(E);
> dc = rph(1+S/2,1+S/2);
>
> srph = ifftshift([A;[B [C;[F dc E];D]]]);
>
> % recreate DFT coeffs, do inverse transform
> nfftimg     = amp .* exp(1i.*srph);
> nimgdata    = ifft2(nfftimg);
>
> % test if its all in real
> sum(real(nimgdata(:)))
> sum(imag(nimgdata(:)))

Subject: randomize phase of Fourier transform and get real data from

From: TideMan

Date: 16 Jan, 2011 19:07:10

Message: 3 of 7

On Jan 13, 9:39 pm, "Diederick C. Niehorster" <dcni...@gmail.com>
wrote:
> Dear all,
>
> I'm puzzling on the following problem:
> I want to randomize the phase spectrum of the Fourier transform of an
> image, but leave the power spectrum untouched, to get a randomized
> image with the same power spectrum as my reference image.
>
> I know that for the inverse transformed to be real, there has to be a
> certain symmetry to the Fourier spectrum, as I use a random phase, I
> need to recreate this. (Is conjugate symmetry the right word?).
>
> I appear to fail however, with a significant part of my image ending
> up in the imaginary part (at least its pretty random, on one run
> almost all is imaginary, at other times its pretty close to all real).
> So I'm probably recreating something wrongly.
>
> Could you tell me where my thinking is in error so I can correct (and
> learn!)?
> Thank you very much for your time,
> Best,
> Dee
>
> % get image from matlab that we all have, see
> %http://blogs.mathworks.com/steve/2006/10/17/the-story-behind-the-matl...
> defimage = pow2(get(0,'DefaultImageCData'),47);
> highbit = 51;
> lowbit = 47;
> numbits = highbit - lowbit + 1;
> b = bitshift(a,-lowbit);
> b = fix(b);
> b = bitand(b,bitcmp(0,numbits));
> theimage = b/max(b(:));
>
> %% here's the start of my code:
> % do fft
> fftimg  = fft2(theimage);
> amp     = abs(fftimg);
> ph      = angle(fftimg);
>
> % random phase
> rph     = 2*pi*rand(size(theimage))-pi;
> % make conjugate symmetric
> %
> % AAAAAAAAAAAAAAAAA
> % B
> % B       C
> % B
> % B   F   dc   E
> % B
> % B       D
> % B
>
> phasemat = fftshift(rph);
>
> S = size(rph,1);
>
> A = rph(1,:);
> B = rph(2:S,1);
> C = rph(2:S/2,2:end);
> D = -flipud(C);
> E = rph(1+S/2,(2+S/2):S);
> F = -fliplr(E);
> dc = rph(1+S/2,1+S/2);
>
> srph = ifftshift([A;[B [C;[F dc E];D]]]);
>
> % recreate DFT coeffs, do inverse transform
> nfftimg     = amp .* exp(1i.*srph);
> nimgdata    = ifft2(nfftimg);
>
> % test if its all in real
> sum(real(nimgdata(:)))
> sum(imag(nimgdata(:)))

Search "phase randomisation" (note spelling) in this newsgroup.

Subject: randomize phase of Fourier transform and get real data from

From: Diederick C. Niehorster

Date: 17 Jan, 2011 06:03:00

Message: 4 of 7

On Jan 17, 3:07 am, TideMan <mul...@gmail.com> wrote:
> On Jan 13, 9:39 pm, "Diederick C. Niehorster" <dcni...@gmail.com>
> wrote:
>
> Search "phase randomisation" (note spelling) in this newsgroup.- Hide quoted text -
>

Hi TideMan,

Thank you, those posts were helpful, showing me a more efficient and
simple way of randomizing the phase spectrum. However, with this
method too, my resulting image is no longer real (having pharand
output its complex result, not just the real part), testing with:
sum(abs(real(nimgdata(:))))
sum(abs(imag(nimgdata(:))))
The energy on the real and the imaginary axes is about equal usually.
This is not something I worry about by now, but something I would like
to understand. How come the output of the inverse Fourier transform is
no longer purely real (barring roundoff error) when the phase is
randomized, even though the symmetry of the Fourier spectrum is
preserved? I hope someone can elp me understand the Fourier transform
better, I'm not able to explain this.

Thank you very much for your help already!
Best,
Dee

Subject: randomize phase of Fourier transform and get real data from

From: TideMan

Date: 17 Jan, 2011 07:12:59

Message: 5 of 7

On Jan 17, 7:03 pm, "Diederick C. Niehorster" <dcni...@gmail.com>
wrote:
> On Jan 17, 3:07 am, TideMan <mul...@gmail.com> wrote:
>
> > On Jan 13, 9:39 pm, "Diederick C. Niehorster" <dcni...@gmail.com>
> > wrote:
>
> > Search "phase randomisation" (note spelling) in this newsgroup.- Hide quoted text -
>
> Hi TideMan,
>
> Thank you, those posts were helpful, showing me a more efficient and
> simple way of randomizing the phase spectrum. However, with this
> method too, my resulting image is no longer real (having pharand
> output its complex result, not just the real part), testing with:
> sum(abs(real(nimgdata(:))))
> sum(abs(imag(nimgdata(:))))
> The energy on the real and the imaginary axes is about equal usually.
> This is not something I worry about by now, but something I would like
> to understand. How come the output of the inverse Fourier transform is
> no longer purely real (barring roundoff error) when the phase is
> randomized, even though the symmetry of the Fourier spectrum is
> preserved? I hope someone can elp me understand the Fourier transform
> better, I'm not able to explain this.
>
> Thank you very much for your help already!
> Best,
> Dee

Without seeing your code, I cannot explain it either.
Perhaps the signal you're inputting to the FFT is complex, not real?
But more likely, you've got a bug in your code.

Subject: randomize phase of Fourier transform and get real data from

From: Diederick C. Niehorster

Date: 21 Jan, 2011 04:12:05

Message: 6 of 7

On Jan 17, 3:12 pm, TideMan <mul...@gmail.com> wrote:
> On Jan 17, 7:03 pm, "Diederick C. Niehorster" <dcni...@gmail.com>
> wrote:
> > However, with this
> > method too, my resulting image is no longer real (having pharand
> > output its complex result, not just the real part), testing with:
> > sum(abs(real(nimgdata(:))))
> > sum(abs(imag(nimgdata(:))))
> > The energy on the real and the imaginary axes is about equal usually.
> > This is not something I worry about by now, but something I would like
> > to understand. How come the output of the inverse Fourier transform is
> > no longer purely real (barring roundoff error) when the phase is
> > randomized, even though the symmetry of the Fourier spectrum is
> > preserved? I hope someone can elp me understand the Fourier transform
> > better, I'm not able to explain this.
>
> > Thank you very much for your help already!
> > Best,
> > Dee
>
> Without seeing your code, I cannot explain it either.
> Perhaps the signal you're inputting to the FFT is complex, not real?
> But more likely, you've got a bug in your code.- Hide quoted text -
>
> - Show quoted text -

Dear Tideman,

Sorry for the slow reply, had a (M.Phil) thesis defense in between.
Thank you for your continued interest.

There was a bug in my code indeed. I started with the pharand function
I found on the forum with your help, but edited it to do a fft2, not
an fft. For the output of fft2 the complex conjugate symmetry doesn't
hold so the pharand method is wrong for that, only figured that out
now. My edit was thus in error.

The reason I am doing a 2D fft is that I want to phase randomize an
image (while preserving the power spectrum). Only doing it in one
direction (with pharand) generates directional noise, while I'd like
it to be randomized in both dimensions. It turns out I'm probably
approaching this fundamentally wrong as the pharand method is not
power spectrum preserving at all either (in my tests at least), if run
over a single dimension.

In my original approach I worked with generating random phases and
recreating their symmetry (following the symmetry in the phases after
fft2), with those phases I reconstructed the Fourier coefficients and
from there simply do ifft2. This preserved the power spectrum, but I
also had some issues with not the whole image being real. The
magnitudes outputted with sum(abs(real(complex_img_data(:)))) and
sum(abs(imag(complex_img_data(:)))) are usually about equal. I think I
am still missing something in how I am creating the phase (maybe I
need the DC parts of the original? I've tried that, made no
difference), or is this a fundamentally wrong approach? Below is the
code that I used to make my random phase spectrum (only works for
square matrices with even sizes).

Thank you very much for your help!
Best,
Dee

% load image and get rid of some hard lines in it
imdata = imread('ngc6543a.jpg'); % comes with matlab
imdata = imdata(1:600,:,1);
imdata(550:end,:) = 0;
imdata(:,590:end) = 0;
figure,imagesc(imdata)

% do fft
fftimg = fft2(imdata);
amp = abs(fftimg);
% get random phase
randphase = rand(size(imdata))*2*pi;
randphase = FourierSymmPhase(randphase);
% reassemble DFT coefficients
nfftimg= amp .* exp(1i.*randphase);
nimdata= ifft2(nfftimg);
% test if image is all real
sum(abs(real(nimdata(:))))
sum(abs(imag(nimdata(:))))

figure,imagesc(real(nimdata))

function symm_mat = FourierSymmPhase(phasemat)
%%
% AAAAAAAAAAAAAAAAA
% B
% B C
% B
% B F dc E
% B
% B D
% B


phasemat = fftshift(phasemat);

S = size(phasemat,1);

A = phasemat(1,:);
B = phasemat(2:S,1);
C = phasemat(2:S/2,2:end);
D = -flipud(C);
E = phasemat(1+S/2,(2+S/2):S);
F = -fliplr(E);
dc = phasemat(1+S/2,1+S/2);

symm_mat = ifftshift([A;[B [C;[F dc E];D]]]);

Subject: randomize phase of Fourier transform and get real data from

From: Ana

Date: 22 Jun, 2013 02:40:28

Message: 7 of 7

"Diederick C. Niehorster" <dcnieho@gmail.com> wrote in message <697e278a-27f8-47c0-9690-ca8399d10a94@21g2000prv.googlegroups.com>...
> On Jan 17, 3:12 pm, TideMan <mul...@gmail.com> wrote:
> > On Jan 17, 7:03 pm, "Diederick C. Niehorster" <dcni...@gmail.com>
> > wrote:
> > > However, with this
> > > method too, my resulting image is no longer real (having pharand
> > > output its complex result, not just the real part), testing with:
> > > sum(abs(real(nimgdata(:))))
> > > sum(abs(imag(nimgdata(:))))
> > > The energy on the real and the imaginary axes is about equal usually.
> > > This is not something I worry about by now, but something I would like
> > > to understand. How come the output of the inverse Fourier transform is
> > > no longer purely real (barring roundoff error) when the phase is
> > > randomized, even though the symmetry of the Fourier spectrum is
> > > preserved? I hope someone can elp me understand the Fourier transform
> > > better, I'm not able to explain this.
> >
> > > Thank you very much for your help already!
> > > Best,
> > > Dee
> >
> > Without seeing your code, I cannot explain it either.
> > Perhaps the signal you're inputting to the FFT is complex, not real?
> > But more likely, you've got a bug in your code.- Hide quoted text -
> >
> > - Show quoted text -
>
> Dear Tideman,
>
> Sorry for the slow reply, had a (M.Phil) thesis defense in between.
> Thank you for your continued interest.
>
> There was a bug in my code indeed. I started with the pharand function
> I found on the forum with your help, but edited it to do a fft2, not
> an fft. For the output of fft2 the complex conjugate symmetry doesn't
> hold so the pharand method is wrong for that, only figured that out
> now. My edit was thus in error.
>
> The reason I am doing a 2D fft is that I want to phase randomize an
> image (while preserving the power spectrum). Only doing it in one
> direction (with pharand) generates directional noise, while I'd like
> it to be randomized in both dimensions. It turns out I'm probably
> approaching this fundamentally wrong as the pharand method is not
> power spectrum preserving at all either (in my tests at least), if run
> over a single dimension.
>
> In my original approach I worked with generating random phases and
> recreating their symmetry (following the symmetry in the phases after
> fft2), with those phases I reconstructed the Fourier coefficients and
> from there simply do ifft2. This preserved the power spectrum, but I
> also had some issues with not the whole image being real. The
> magnitudes outputted with sum(abs(real(complex_img_data(:)))) and
> sum(abs(imag(complex_img_data(:)))) are usually about equal. I think I
> am still missing something in how I am creating the phase (maybe I
> need the DC parts of the original? I've tried that, made no
> difference), or is this a fundamentally wrong approach? Below is the
> code that I used to make my random phase spectrum (only works for
> square matrices with even sizes).
>
> Thank you very much for your help!
> Best,
> Dee
>
> % load image and get rid of some hard lines in it
> imdata = imread('ngc6543a.jpg'); % comes with matlab
> imdata = imdata(1:600,:,1);
> imdata(550:end,:) = 0;
> imdata(:,590:end) = 0;
> figure,imagesc(imdata)
>
> % do fft
> fftimg = fft2(imdata);
> amp = abs(fftimg);
> % get random phase
> randphase = rand(size(imdata))*2*pi;
> randphase = FourierSymmPhase(randphase);
> % reassemble DFT coefficients
> nfftimg= amp .* exp(1i.*randphase);
> nimdata= ifft2(nfftimg);
> % test if image is all real
> sum(abs(real(nimdata(:))))
> sum(abs(imag(nimdata(:))))
>
> figure,imagesc(real(nimdata))
>
> function symm_mat = FourierSymmPhase(phasemat)
> %%
> % AAAAAAAAAAAAAAAAA
> % B
> % B C
> % B
> % B F dc E
> % B
> % B D
> % B
>
>
> phasemat = fftshift(phasemat);
>
> S = size(phasemat,1);
>
> A = phasemat(1,:);
> B = phasemat(2:S,1);
> C = phasemat(2:S/2,2:end);
> D = -flipud(C);
> E = phasemat(1+S/2,(2+S/2):S);
> F = -fliplr(E);
> dc = phasemat(1+S/2,1+S/2);
>
> symm_mat = ifftshift([A;[B [C;[F dc E];D]]]);

You can get rid of the negligible imaginary part by using real... (http://blogs.mathworks.com/steve/2010/07/16/complex-surprises-from-fft/)

Tags for this Thread

No tags are associated with this thread.

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us