I cannot match the calculated from formula FFT with the MATLAB answer though they should be the same

I have a formula for a Soliton. I have a formula for its FFT but I cannot get the calculated version to match the theoretical formula. It appears that Engineers, Mathematicians, Physicists and DSP all normalise in different ways!! I need to know how to correct my code to get them synchronised. The result however is a massive overstatement by the theoretical value. If nothing else I do not have a grip on the normalisation. An explanation of how I am misusing/misunderstanding the FFT process would be much appreciated.

%making a dummy soliton
            close all;
            a = 1; g = 1; v = 0;% a - amplitude, g - gamma, v - velocity.
            %b = sqrt(g/2)*a;
            xR = 10; %xRange
            %s1 = @(x) a*exp(1i*(v/2)*x).*sech(b*x); %Soliton
            %f1 = @(w) 1/sqrt(pi)*a*sech(pi*(w+v/2)/(sqrt(2)*a))/abs(a);
            s1 = @(x) a*sech(a/sqrt(2)*x); %Soliton
            f1 = @(w) sqrt(pi)*a*sech(pi*w/(sqrt(2)*a))/abs(a);
            h = 0.01; %dx
            x = -xR:h:xR;
            S = s1(x);
            % figure(1);
            % plot(x,abs(S)); grid on;
            figure(2);
            NFFT = 2^nextpow2(length(x));
            X = fftshift(fft(S,NFFT));
            fVals = (-NFFT/2:NFFT/2-1)/NFFT;
            fs = 1/h; fVals = fs*fVals;
            L = length(X);
            Px = (X.*conj(X))/(NFFT*L);
            Am = sqrt(Px);
            plot(fVals,Px,'b',DisplayName='Power Spectral Density');grid on; hold on;
            plot(fVals,Am,'r', DisplayName='Magnitude')
            xlabel('Frequency');
            ylabel('Magnitude');
            xlim([-5 5])
            legend;
            F = f1(x);
            plot(x,F,'k',DisplayName='Theoretical Amplitude'); grid on;

Answer 1

[Paul]

Hi Jonathan

Let's define the input signal symbolically

clear
syms a x w real
s1(x) = a*sech(a/sqrt(2)*x)

s1(x) = 

Its Continous Time Fourier Transform (CTFT) is with frequency variable w in rad/sec

f1(w) = simplify(fourier(s1(x),x,w))

f1(w) = 

This doesn't match the equation in the code in the Question. There is a difference of sqrt(2) in the numerator, and the sign will be different if a < 0.

sqrt(sym(pi))*a*sech(sym(pi)*w/(sqrt(2)*a))/abs(a)

ans = 

Use the specific value of a

a = 1;
s1(x) = subs(s1(x))

s1(x) = 

f1(w) = subs(f1(w))

f1(w) = 

Convert to anonymous functions for numerical computation

s1 = matlabFunction(s1)

s1 = function_handle with value:
    @(x)1.0./cosh((sqrt(2.0).*x)./2.0)

f1 = matlabFunction(f1)

f1 = function_handle with value:
    @(w)(sqrt(2.0).*pi)./cosh((sqrt(2.0).*w.*pi)./2.0)

Define the evaluation values for x and plot

h = 0.01; %dx
xR = 10; %xRange
x = -xR:h:xR;
S = s1(x);
figure
plot(x,S),grid

The range of x covers the significant nonzero portion of s1.

Compute the DFT. Because s1 is evaluated for x < 0, we would need to shift the input before taking the DFT if we want the phase of the DFT to match the phase of the CTFT, which is zero. But if we only care about amplitude we don't need to bother and can proceed as

NFFT = 2^nextpow2(length(x));
X = fftshift(fft(S,NFFT));

I'll rename to wVals because we are working in rad/sec. We could work in Hz with appropriate substitution of f for w in f1.

%fVals = (-NFFT/2:NFFT/2-1)/NFFT;
%fs = 1/h; fVals = fs*fVals;
wVals = (-NFFT/2:NFFT/2-1)/NFFT*2*pi/h;

Plot the amplitude of the CTFT

figure
plot(wVals,abs(f1(wVals)))
hold on

Plot the ampltude of the DFT. With the default definitions in Matlab, the DFT has to be scaled by h to approximate the CTFT.

stem(wVals,abs(X)*h)
xlim([-10 10])
xlabel('rad/sec')

[Jonathan Wharrier]

What can I say... I think I understood this AND it answers the question. Truly grateful!

Jo