Why Does ifourier Not Return Equivalent Results for Equivalent Transforms?

Consider two expressions of a Fourier transform

syms t omega real
HU(omega) = -(exp(-omega*1i)*(exp(omega*1i) - 1)*(omega + 1i))/(omega*(omega^2 + 1))

HU(omega) = 

Y(omega) = (exp(-omega*1i)*1i - 1i)/(omega*(1 + omega*1i))

Y(omega) = 

The transforms are equivalent

simplify(HU(omega)-Y(omega))

ans = 

0

so I would expect the inverse transforms to be equivalent.

Inverse transform of Y(omega)

y(t) = ifourier(Y(omega),omega,t)

y(t) = 

Inverse transform of HU(omega)

hu(t) = ifourier(HU(omega),omega,t)

hu(t) = 

hu(t) looks complicated, but the primary concern is that Dirac delta in the first term in the numerator.

Simplify both expressions

s = @(z) simplify(rewrite(simplify(z),'heaviside'));
y(t) = s(y(t))

y(t) = 

hu(t) = s(hu(t))

hu(t) = 

Sure enough, both time domain signals are the same with the exception that hu(t) includes that Dirac delta term.

Show that they are the same, except for the Dirac delta, by subtracting and using fplot that ignores Dirac delta terms:

figure
fplot(y(t)-hu(t),[-10,10]),ylim([-1e-16,1e-16])

Is this a bug, or is there a problem in the code, or am I misunderstanding what the results should be?

[Walter Roberson]

Taking a round trip fourier(ifourier(X)) is apparently not producing the same answer for hu -- they differ by a constant of 1

syms t omega real
HU(omega) = -(exp(-omega*1i)*(exp(omega*1i) - 1)*(omega + 1i))/(omega*(omega^2 + 1))

HU(omega) = 

Y(omega) = (exp(-omega*1i)*1i - 1i)/(omega*(1 + omega*1i))

Y(omega) = 

figure(1)
y(t) = ifourier(Y(omega),omega,t)

y(t) = 

hu(t) = ifourier(HU(omega),omega,t)

hu(t) = 

yy = rewrite(simplify(y), 'heaviside')

yy(t) = 

symvar(yy)

ans = 

t

YY = fourier(yy)

YY = 

YY = subs(YY, symvar(YY), omega)

YY = 

YY - Y

ans(omega) = 

simplify(ans)

ans(omega) = 

0

simplify(expand(ans))

ans(omega) = 

0

C1 = collect(ans, omega)

C1(omega) = 

0

symvar(C1)

ans = 

ω

fplot(C1, [0 5])

figure(2)
huhu = rewrite(simplify(hu), 'heaviside')

huhu(t) = 

HUHU = fourier(huhu)

HUHU = 

HUHU = subs(HUHU, symvar(HUHU), omega)

HUHU = 

HUHU - HU

ans(omega) = 

simplify(ans)

ans(omega) = 

simplify(expand(ans))

ans(omega) = 

C2 = collect(ans, omega)

C2(omega) = 

symvar(C2)

ans = 

ω

fplot(C2, [0 5])

[Paul]

Yes, the Fourier transform of the Dirac delta is 1, so the -1 resulting from HUHU - HU is the Fourier transform of the -dirac(t) that popped up in hu(t) and stuck around in huhu(t). Looks like a bug to me ....

FWIW, I believe that y(t) is correct and hu(t) is incorrect.

[David Goodmanson]

Hi Paul,

That would appear to be the case, especially since (putting things slightly differently and ignoring multiplicative constants),

  lim{w-->inf} hu(w)
= lim{w-->inf} Ftrans(hu(t))     % supposed hu(t)
= lim{w-->inf} Ftrans(delta(t) + sum of heavisides)
= lim{w-->inf} (1 + O(1/w)) = 1                      *

whereas

hu(w) = y(w) is O(1/w^2) as w-->inf 

* Using O(1/w) for one heaviside and still O(1/w) for a sum of heavisides is a simple upper bound to illustrate the point. For a particular sum of heavisides the actual behavior could fall off faster, as it does here.

[Paul]

I've opened a support case. I'll report back here when it gets resolved.

Answer 1

[Paul]

Tech Support confirmed a bug (and it's quite pernicious IMO).

Their suggested workaround for this particular problem is to simplify HU(omega), but it needs more than the default number of steps. In this case, 200 seems to do the trick (I did not try any other number of steps)

syms t omega real
HU(omega) = -(exp(-omega*1i)*(exp(omega*1i) - 1)*(omega + 1i))/(omega*(omega^2 + 1))

HU(omega) = 

simplify(HU(omega))

ans = 

HU(omega) = simplify(HU(omega),200)

HU(omega) = 

hu(t) = simplify(rewrite(ifourier(HU(omega),omega,t),'heaviside'))

hu(t) = 

My experience with fourier is that it prefers expressions in exp over sin and/or cos and I would have assumed that ifourier prefers the same. The irony is that if I had been given HU(omega) as shown above I probably would have rewritten it in terms of exp as

simplify(rewrite(HU(omega),'exp'))

ans = 

and then suffered from the wrong result after taking the inverse.

Tech Support said the bug is related to the fact that HU(omega) has common factor in its numerator and denominator, which is (1 - 1i*omega), though I have no idea why that would cause a problem.

[num,den] = numden(HU(omega));
term = (1 - 1i*omega);
num = simplify(num/term);
den = simplify(den/term);
HU(omega) = num/den

HU(omega) = 

hu(t) = simplify(rewrite(ifourier(HU(omega),omega,t),'heaviside'))

hu(t) = 

Of course, how one would know that's a problem for a somewhat complicated expression is unclear.

Also, a common factor doesn't appear to be a problem in general. For example

Z(omega) = 1/(1 + 1i*omega);
z(t) = ifourier(Z(omega),omega,t)

z(t) = 

Define Z(omega) with a common term in the numerator and denominator

Z(omega) = (1 - 1i*omega)/expand((1 + 1i*omega)*(1 - 1i*omega))

Z(omega) = 

z(t) = simplify(ifourier(Z(omega),omega,t))

z(t) = 

Same answer.