清空筛选条件
清空筛选条件

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

8 次查看（过去 30 天）
Jonathan Wharrier
Jonathan Wharrier 2023-3-29
编辑： Paul 2023-3-30
Ran in:
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;

请先登录，再回答此问题。

采纳的回答

Paul
Paul 2023-3-30
编辑：Paul 2023-3-30
Ran in:
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')

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Discrete Fourier and Cosine Transforms 的更多信息

标签

产品


版本

R2023a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by