Failing to plot a fourier transform

Hello,

I've been trying to understand the functioning of modulation and demodulation using a cosine wave, and study how it works in time and frequency domains. Thus, I've put into matlab what the book does, to see it visually.

 syms x;
>> y1 = power(sinc(x),2);
>> y2 = cos(2*pi*x);
>> z1 = times(y1,y2)
>> f1 = fourier(y1);
>> f2 = fourier(y2);
>> f3 = fourier(y3);

However, when using fplot to plot the first three (so, the sinc^2, the cosine wave and then the modulated signal z1 = sinc^2(x) * cos(2πt)) everything goes smoothly. Wanting to see whether or not this time-domain function does indeed equal the frequency-domain plots the book shows, I transformed every signal. In theory, f1 should be a Tri, f2 should be two dirac impulses, and f3 should be the convolution of f1 and f2. Here's the issue: when using fplot to plot the frequency-domain functions, I get the following error:

Error using fplot>singleFplot
Input must be a function or functions of a single variable.

Error in fplot>@(f1,f2)singleFplot(cax,{f1,f2},limits,extraOpts,args) (line 208)
        hObj = cellfun(@(f1,f2) singleFplot(cax,{f1,f2},limits,extraOpts,args),fn{1},fn{2},'UniformOutput',false);

Error in fplot>vectorizeFplot (line 208)
        hObj = cellfun(@(f1,f2) singleFplot(cax,{f1,f2},limits,extraOpts,args),fn{1},fn{2},'UniformOutput',false);

Error in fplot (line 166)
    hObj = vectorizeFplot(cax,fn,limits,extraOpts,args);

So, I tried to integrate over x (which, in this example, is time):

>> f1 = fourier(y1,x);
>> f2 = fourier(y2,x);
>> f3 = fourier(z1,x);

Which then led to the following error:

Warning: Error updating ParameterizedFunctionLine.

 The following error was reported evaluating the function in FunctionLine update: Unable to convert expression containing
 remaining symbolic function calls into double array. Argument must be expression that evaluates to number.

Is it possible to plot the frequency-domain functions? If so, what am I doing wrong?

I apologise if the code looks very simple and lackluster but this is my first time using matlab.

Thanks in advance

[Walter Roberson]

 syms x;
 y1 = power(sinc(x),2)

y1 = 

 y2 = cos(2*pi*x)

y2 = 

 z1 = times(y1,y2)

z1 = 

 f1 = fourier(y1)

f1 = 

 f2 = fourier(y2)

f2 = 

 fplot(z1, [0 2*pi])

 %dirac deltas are either 0 or infinite so no useful plot
 fplot(f2, [0 100])

 fplot(f1, [0 100]) 

Warning: Error updating FunctionLine.

The following error was reported evaluating the function in FunctionLine update: Unable to convert expression containing remaining symbolic function calls into double array. Argument must be expression that evaluates to number.

That is caused by the fact that fourier() does not know how to transform cos(x)/x^2 . Wolfram Alpha does know how to transform sinc(x)^2, getting sqrt(pi)/2 times a triangle distribution, or an expression involving five heaviside functions.

[Lorenzo Romagnoli]

So, if I am understanding this correctly, I cannot plot the frequency-domain functions in matlab, and thus I need to rely on WolframAlpha?

[Walter Roberson]

You could switch from symbolic fourier transform to discrete fourier transform fft

[Lorenzo Romagnoli]

What changes if I do that?

[Walter Roberson]

Choose particular x to evaluate the function at; use regular intervals. fft() the function values. This will give you amplitudes at frequencies determined by the spacing in your x values.

[Walter Roberson]

Hmmm, except that fft assumes that the signal is periodic. That can make a significant difference in the plot.

Answer 1

[Paul]

Hi Lorenzo,

On the face of it, the Symbolic Math Toolbox treatment of (normalized) sinc(x) is disappointing. Perhaps there are some nuances I don't understand that makes it more difficult to handle than I imagine it to be.

Evaluate sinc(x) at the origin

sinc(sym(0))

ans = 

1

So for, so good. Define a symfun in terms of sinc(x)

syms x real
y0(x) = sinc(x)

y0(x) = 

For some reason, the SMT replaces sinc(x) with the expression shown, instead of just carrying sinc along, like cos or exp. Howvever, y0(x) is not actually equal to sinc(x) because x = 0 has to be excluded from the domain of y0(x)

try y0(0)
catch ME
    disp('error')
end

error

Clearly, y0(x) is not equal to sinc(x) because sinc(x) is well-defined at the origin.

However, for some reason, isAlways evaluates to true

isAlways(y0(x) == sinc(x))

ans = logical
   1

which, according to isAlways, should only be the case "if cond holds true for all possible values of the symbolic variables in cond including all assumptions on the variables" (emphasis added). This result is a bug, or the documentation is incorrect.

Taking the FT of y0

syms w real
rect(w) = rectangularPulse(w);
f0(w) = fourier(y0(x),x,w)

f0(w) = 

which is the same as rect(w/(2pi)), as expected (at least for plotting, I wasn't able to do the symbolic manipulations to show it)

fplot(f0(w))

fplot([f0(w) - rect(w/2/sym(pi))])

Moving on to sinc(x)^2 ....

As you saw, the SMT can't find the FT of sinc(x)^2

y1(x) = sinc(x)^2

y1(x) = 

f1(x) = fourier(y1(x),x,w)

f1(x) = 

It seems like this would be much easier if sinc was carried along on the rhs and then fourier would have the FT of sinc and sinc^2 stored in a table, like I imagine it stores many others. But even though it doesn't, it would seem that fourier() could still get the correct result using the "multiply by sin" property of the FT (again, I'm assuming that part of the magic of fourier is that it uses properties of the FT to derive results for general expressions).

f1(w) = fourier(1/x^2/sym(pi)^2,x,w);              % FT of 1/(x*pi)^2
f1(w) = (f1(w - sym(pi)) - f1(w + sym(pi)))/2/1i;  % FT after multiply by sin(pi*x)
f1(w) = (f1(w - sym(pi)) - f1(w + sym(pi)))/2/1i;  % FT after multiply by sin(pi*x) again
f1(w)

ans = 

This result is exactly what we expect, i.e., a triangular pulse

tri(w) = triangularPulse(w);
fplot(f1(w),[-8 8])

fplot(f1(w) - tri(w/2/sym(pi)),[-8 8])

As an aside, the ifourier of tri(x/(2pi)) is

g1(x) = ifourier(tri(w/2/sym(pi)),w,x)

g1(x) = 

which is an expression for sinc(x)^2 (except at the origin)

simplify(g1(x)-sinc(x)^2)

ans = 

0

This result suggests another way to compute f1

f1(w) = simplify(fourier(rewrite(y1(x),'exp'),x,w),100)

f1(w) = 

which is the same as above (don't know why simplify can't do better).

Moving on to y3 (which was z3 in the Question), it's not suprising given the above results that fourier fails. However, starting with f1(w) we could derive f3(w) using the "multiply by cos" property of the FT, or maybe the rewrite in terms of exp would do the trick.

As for y2(x), fourier does come up with the correct result

y2(x) = cos(2*sym(pi)*x);
f2(w) = fourier(y2(x),x,w)

f2(w) = 

As previously noted by @Walter Roberson, fplot ignores dirac, so the fplot of f2 shows as zero everywhere. I think it may be feasible for a general expression to strip out the dirac functions, fplot the leftover expression, and then use stem to add vertical arrows to the plot for the diracs.

For y0, y1, and y3, we could use discrete-time techniques to approximate their continuous-time FTs. We could do the same for y2 if we a) carefully select the sampling period and b) carefully select the number of samples, and c) understand the relationship between discrete- and continuous-time Fourier series coefficients and the continuous-time Fourier transform (I'm using "time" here in a generic sense of the independent variable, which in this case is x).

Summary: It seems like the SMT treatment of sinc(x) could be improved.