0
$\begingroup$

I'm trying to implement the code in this tutorial in MATLAB using MATLAB's ifft2() function instead of the algorithm used in the shader (looks to be some type of Radix-2 algorithm). It is used to compute the surface of an ocean using compute shaders in Unreal Engine.

I've confirmed I'm computing FFTGrid.DZ correctly, but when I take the inverse fast Fourier transform of it (D = ifft2(hTilde)) I get a lot of noise in the corners as shown below and the results don't look correct. Are there common causes for this?

enter image description here

Thanks!

EDIT

The MATLAB 2023a code I'm executing is:

clc; clear all; close all;

%%-- FX_OceanWater_SetInitials --%%
GridSize = 64;
HalfGridSize = GridSize/2;

% PerCascadeParameters
WindDirectionality = 1;
WindDirectionality = 1-WindDirectionality;
A = 240;
PatchLength = 2000; % [m]
Chop = 1.5;
ShortWaveCutoff = 30;
LongWaveCutoff  = 0.04;
WindTighten = 1;

% Misc
g = 9.8; % [m/s^2]
RepeatPeriod = 1000;
Time = 0.0125;
AnimationTime = mod(Time,RepeatPeriod);
BaseFreq = 2*pi/RepeatPeriod;

% WindControl
windAngle = 90;
WindDirection = [sind(windAngle);cosd(windAngle)];
WindSpeed = 44;


%%-- FX_OceanWater_SetGrids--%%
SpectrumGrid    = cell(GridSize,GridSize);
FFTGrid.DX_DY   = cell(GridSize,GridSize);
FFTGrid.DZ      = cell(GridSize,GridSize);


%%-- FX_OceanWater_PopulateSpectrum --%%
threadX  = 0:1:GridSize-1;
threadY = 0:1:GridSize-1;

% Complex Amplitudes
rng(1);
randNums = ones(4,1);%rand(4,1)*2*pi;
ComplexAmplitude = [randNums(1)*sin(randNums(1));randNums(2)*cos(randNums(2));randNums(3)*sin(randNums(3));randNums(4)*cos(randNums(4))];

% HLSL
for ii = 1:length(threadX)
    for jj = 1:length(threadY)
        % Retrive WaveVector from thread index.
        Position = [threadX(ii);threadY(jj)];
        WaveVector = [Position(1)-HalfGridSize; Position(2)-HalfGridSize];
        WaveVector = WaveVector * 2 * pi;
        
        % Our parameters for each of cascade are stored as components of a 4D vector.
        % We access them, treating 4D vector as array of 4 scalars,
        % And our Z component of thread index is the index of cascade, current compute thread belongs to.
        WaveVector = WaveVector/PatchLength;
                
        % Calculate magnitude of WaveVector
        k = norm(WaveVector,2);

        % Result is magnitude of positive and negative spectrums
        Result = [0;0];

        % Handle zero length case
        if k > 0.0001
            % Calculate direction of wave vector
            k_norm = WaveVector/k;
            
            % Calculate dampening factor for wave directions not aligning with the wind
            % Positive spectrum is in R component, negative spectrum is in G component.
            WindFactor(1) = dot( k_norm, WindDirection);
            WindFactor(2) = dot(-k_norm, WindDirection);
            WindFactorAbs = abs(WindFactor);
            WindFactorAbs = WindFactorAbs.^WindTighten;

            % Reduce magnitude of the waves, travelling in negative direction
            if WindFactor(1) <= 0
                WindFactorAbs(1) = WindFactorAbs(1)*WindDirectionality;
            end
            if WindFactor(2) <= 0
                WindFactorAbs(2) = WindFactorAbs(2)*WindDirectionality;
            end
            WindFactor = WindFactorAbs;

            % Phillips Ocean spectrum calculation
            L = WindSpeed * WindSpeed / g;
            UpperPart = exp(-1.0 / ((k * L) * (k * L)));
            Spectrum = A * UpperPart / (k * k * k * k);

            % Dampen waves, shorter than user controlled threshold
            Spectrum = Spectrum * exp(-(k*k) * ShortWaveCutoff);
            if k < LongWaveCutoff
                Spectrum = 0;
            end

            Result = sqrt(Spectrum * WindFactor);
            Result = Result / sqrt(2);
        end

        % Lastly, obtain complex amplitude/phases.
        PositiveAndNegativeSpectrum = [Result(1); Result(1); Result(2); Result(2)];
        PositiveAndNegativeSpectrum = PositiveAndNegativeSpectrum.*ComplexAmplitude;

        % And obtain conjugate of negative spectrum by flipping its imaginary part.
        PositiveAndNegativeSpectrum(4) = -1*PositiveAndNegativeSpectrum(4);

        % Prepare index of the grid, at which result will be stored and write the result.
        SpectrumGrid{ii,jj} = PositiveAndNegativeSpectrum;
    end
end


%%-- FX_OceanWater_TimeStep --%%
% HLSL
for ii = 1:length(threadX)
    for jj = 1:length(threadY)
        % Retrive WaveVector from thread index
        WaveVector = [threadX(ii)-HalfGridSize;threadY(jj)-HalfGridSize];
        WaveVector = WaveVector * 2 * pi;
        WaveVector = WaveVector/PatchLength;
        
        % Calculate magnitude of WaveVector
        k = norm(WaveVector,2);
        
        % Dispersion. Angular frequency from wave vector
        Freq = sqrt(g * k);
        
        % Quantize frequency to multiple of base frequency
        Freq = floor(Freq / BaseFreq) * BaseFreq;
        Phase = Freq * AnimationTime;

        % Load initial spectrum
        h0 = SpectrumGrid{ii,jj};

        % Calculate exponents for positive and negative spectrums
        fourier_amp         = h0(1:2);
        fourier_amp_conj    = h0(3:4);
        SineCosine = [sin(Phase);cos(Phase)];

        exponent = [SineCosine(2); SineCosine(1)];
        exponent_inv = [SineCosine(2); -SineCosine(1)];

        % Resulting displacements
        DispZ=[0;0];
        DispX=[0;0];
        DispY=[0;0];

        % Handle zero length case.
        if k > 0.000001
            %Complex multiplication of positive spectrum by exponent
            c0 = jMul(fourier_amp,exponent);

            %Complex multiplication of negative spectrum by inverse exponent
            c1 = jMul(fourier_amp_conj,exponent_inv);

            %Complex addition of positive and negative parts
            DispZ = c0 + c1;

            % Calculate horizontal displacements by projecting vertical displacement on components of wave direction
            dx = [0;WaveVector(1)/k]*Chop;
            DispX = jMul(DispZ,dx);
            
            dy = [0;  WaveVector(2)/k]*Chop;
            DispY = jMul(DispZ,dy);
        end
        FFTGrid.DX_DY{ii,jj} = [DispX;DispY];
        FFTGrid.DZ{ii,jj} = [DispZ];
    end
end



%% Compute Z
% Convert FFTGrid to complex array
for ii = 1:length(threadX)
    for jj = 1:length(threadY)
        hTilde(ii,jj) = FFTGrid.DZ{ii,jj}(1) + ii*FFTGrid.DZ{ii,jj}(2);
    end
end

%- Definitions -%
% -N/2 <= n < N/2
n = -GridSize/2:1:(GridSize/2-1);

% -M/2 <= m < M/2
m = -GridSize/2:1:(GridSize/2-1);

% Discrete spatial points
xPoints = n*PatchLength/GridSize;
zPoints = m*PatchLength/GridSize;

% x: Spatial Points
[xGrid,zGrid] = meshgrid(xPoints,zPoints);

% Compute iFFT
D = ifft2(hTilde);

%-- Draw Surface --%
s=surf(xGrid,zGrid,real(D));
s.EdgeColor = 'none';


%% User Functions
function sgn = signGrid(N,M)
    [x,y] = meshgrid(1:M,1:N) ;
    sgn = ones( M, N ) ;
    sgn(mod(x+y,2)==0) = -1 ;
end

function c = jMul(c0, c1)
    c       = zeros(2,1);
    c(1)    = c0(1)*c1(1) - c0(2)*c1(2);
    c(2)    = c0(1)*c1(2) + c0(2)*c1(1);
end

where hTilde are the Fourier amplitudes of the wave field realization.

Update I believe fftshift does shift the results to the proper location, but it doesn't eliminate or explain the large spikes

surf(xGrid,zGrid,fftshift(real(D))): enter image description here

$\endgroup$
7
  • 3
    $\begingroup$ Are you sure you aren't missing an fftshift? $\endgroup$
    – Baddioes
    Feb 22 at 19:33
  • $\begingroup$ I could be? When are you supposed to use that? $\endgroup$
    – eball
    Feb 22 at 20:23
  • 2
    $\begingroup$ ifft assumes your spectrum is from $[0,2\pi)$. If you generate a spectrum that's from $[-\pi,\pi)$ you need to ifftshift() it prior to performing an IFFT. Ie, if you have a sinusoid with digital frequency $0.1\pi$ and create the spectrum from $[-\pi,\pi)$, it will appear to the ifft that the frequency is actually at $0.9\pi$ unless you ifftshift. $\endgroup$
    – Baddioes
    Feb 22 at 21:20
  • $\begingroup$ Please read about MRE, then edit your post to include one. If we don’t know what FFTGrid.DZ.x is, we don’t know what you’re computing the inverse transform of. $\endgroup$ Feb 23 at 15:24
  • $\begingroup$ I'm not sure how to make the code shorter while being reproducible. I'll include the code I'm running and think about a way to make it shorter $\endgroup$
    – eball
    Feb 23 at 17:46

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.