Failing to plot a fourier transform
6 次查看(过去 30 天)
显示 更早的评论
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
6 个评论
Walter Roberson
2022-11-6
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
2022-11-6
Hmmm, except that fft assumes that the signal is periodic. That can make a significant difference in the plot.
采纳的回答
Paul
2022-11-13
编辑:Paul
2022-12-14
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))
So for, so good. Define a symfun in terms of sinc(x)
syms x real
y0(x) = sinc(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
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))
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)
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
f1(x) = fourier(y1(x),x,w)
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)
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])
g1(x) = ifourier(tri(w/2/sym(pi)),w,x)
which is an expression for sinc(x)^2 (except at the origin)
simplify(g1(x)-sinc(x)^2)
This result suggests another way to compute f1
f1(w) = simplify(fourier(rewrite(y1(x),'exp'),x,w),100)
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)
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.
0 个评论
更多回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Calculus 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!