Hello everyone, I like to ask how do i plot a function that integrates a certain function
Cn = integral(function of the sawtooth * exp(-1i*2*pi*n*f*t), 0, T);
so far this is my working code which only plots the sawtooth wave function. I have made multiple changes to it trying to figure out how to plot Cn. I'm still fairly new to matlab. maybe I'm missing something.
close all; figure; t=linspace(0,3,1500); m=linspace(0,3,1500); T=1; f=1/T; amp=1; sawf=amp/2; k=plot(NaN,NaN); 1i; for n=1:1:500 sawf = sawf - (amp/pi)*(1/n)*(sin(2*pi*n*f*t)); set(k, 'XData',t,'YData',sawf); pause(0.01); end
%cn1= integral(sawf*((cos(2*pi*1*f*t))-(-(1i)*sin(2*pi*1*f*t))),0,T); %cn1= sawf*(sin(2*pi*n*f*t)/(2*pi*n*f))+((i)cos(2*pi*n*f*t)/(2*pi*n*f)); %cn = (1./T)*cn1; %set(k, 'XData',m,'YData',cn);

댓글 수: 3

here's the code in proper format. for some reason cloudflare won't let me post on my pc so i was forced to post using my phone
close all;
figure;
t=linspace(0,3,1500);
m=linspace(0,3,1500);
T=1;
f=1/T;
amp=1;
sawf=amp/2;
k=plot(NaN,NaN);
1i;
for n=1:1:500
sawf = sawf - (amp/pi)*(1/n)*(sin(2*pi*n*f*t));
set(k, 'XData',t,'YData',sawf);
pause(0.01);
end
%cn1= integral(sawf*((cos(2*pi*1*f*t))-(-(1i)*sin(2*pi*1*f*t))),0,T);
%cn1= sawf*(sin(2*pi*n*f*t)/(2*pi*n*f))+((i)cos(2*pi*n*f*t)/(2*pi*n*f));
%cn = (1./T)*cn1;
%set(k, 'XData',m,'YData',cn);
Dyuman Joshi
Dyuman Joshi 2023년 9월 9일
편집: Dyuman Joshi 2023년 9월 9일
integral requires a function handle as an integrand. You have not supplied any function, just a numerical value.
Also, I don't understand how this -
%1
cn1= integral(sawf*((cos(2*pi*1*f*t))-(-(1i)*sin(2*pi*1*f*t))),0,T);
results in this -
%2
cn1= sawf*(sin(2*pi*n*f*t)/(2*pi*n*f))+((i)cos(2*pi*n*f*t)/(2*pi*n*f));
or how did you arrive at line 2.
or was it supposed to be n instead of 1 (next to pi inside the trignometric terms) in line 1?
jvfa
jvfa 2023년 9월 9일
Hi, can you elaborate where i where i got wrong on what you said about integral? because it's really a big part of understanding the matlab function of integration that I lack.
those comments don't worry about them they were just my attempts at understanding the function. %1 was originally n then i tried substituting 1 for n to see what the value is sort of doing trial and error. for %2 i manually integrated the thing by hand, given that i have the value for sawf so i can move it out of the integral sign then just integrated exp(-2i*pi*n*f*t)dt instead

댓글을 달려면 로그인하십시오.

 채택된 답변

Bruno Luong
Bruno Luong 2023년 9월 9일
편집: Bruno Luong 2023년 9월 9일

0 개 추천

The coefficients Cn = 1i/(2*pi*n).
I'm not very good in symbolic calculation, I never own the tollbox license.
syms t
syms n integer
sawtooth = t;
Cn = simplify(int(sawtooth*exp(-1i*2*pi*n*t), t, 0, 1))
Cn = 
C0 = simplify(int(sawtooth, t, 0, 1)) % funny the above general formula is wrong for n=0
C0 = 
Cnfun = matlabFunction(Cn);
n = -10:10;
Cnval = Cnfun(n);
Cnval(n==0) = subs(C0); % no idea why I need to do this
figure
hold on
plot(n, real(Cnval))
plot(n, imag(Cnval))
legend('real','imag')
xlabel('n')
ylabel('C(n)')
t = linspace(0,3,1000)';
f = @(t) mod(t,1)
f = function_handle with value:
@(t)mod(t,1)
fFourier = sum(Cnval.*exp(1i*2*pi*n.*t),2);
figure
plot(t, f(t), 'r', t, fFourier, 'b')
Warning: Imaginary parts of complex X and/or Y arguments ignored.
xlabel('t')
legend('f', 'Fourier approximation')

댓글 수: 5

jvfa
jvfa 2023년 9월 9일
I didn't really understand what you said there, also I tried to run your code it's a fourier approximation plot? I'm trying to plot something else i think
Bruno Luong
Bruno Luong 2023년 9월 9일
"I have a homework that asks me to plot Cn with respect to n"
that's the first plot, Cn vs n
Paul
Paul 2023년 9월 9일
편집: Paul 2023년 9월 9일
Hi Bruno,
Regarding the concern about the general formula being wrong ...
syms t
syms n integer
sawtooth = t;
Don't use simplify to see what's happening (always safer to use sym(pi) )
Cn = int(sawtooth*exp(-1i*2*sym(pi)*n*t), t, 0, 1)
Cn = 
int comes up with a general solution that doesn't take into account the special case of n = 0. However, note that it's sort of the right answer for n = 0 in the sense that
limit(Cn,n,0)
ans = 
This situation comes up all of the time with Fourier integrals. I'm not sure if it's reasonable to expect that int identify all of the special cases for all parameters in the integrand and then provide a result in terms of piecewise that captures each special case. I've often wondered about this issue.
I don't know the inner workings of simplify, but I'll note that we can do this
[num,den] = numden(Cn)
num = 
den = 
num = simplify(num)
num = 
And then
num/den
ans = 
which is standard for the Symbolic toolbox, i.e., n/n will always be replaced with 1, even if n = 0 is a possibility.
C0 = simplify(int(sawtooth, t, 0, 1)) % funny the above general formula is wrong for n=0
C0 = 
What I do is identify the roots of the denominator of the output of int, before doing any simplification, and then do a piecewise and limit. The roots can be found using solve, for example, but here it's easy to do by hand
Cn = piecewise(n == 0, limit(Cn,n,0), Cn)
Cn = 
Then simplify at the end
Cn(n) = simplify(Cn)
Cn(n) = 
The drawback to this approach is that piecewise is not supported by matlabFunction
try
Cnfun = matlabFunction(Cn);
catch ME
ME.message
end
ans = 'Unable to generate code for piecewise for use in anonymous functions. Use option 'File' to write code for piecewise to a file instead.'
So I just evaluate Cn symbolically and convert to double (if numeric is preferred)
n = -10:10;
%Cnval = Cnfun(n);
Cnval = double(Cn(n));
% subs wasn't needed here
% Cnval(n==0) = subs(C0); % no idea why I need to do this
figure
hold on
stem(n, real(Cnval))
stem(n, imag(Cnval))
legend('real','imag')
xlabel('n')
ylabel('C(n)')
t = linspace(0,3,1000)';
f = @(t) mod(t,1);
fFourier = sum(Cnval.*exp(1i*2*pi*n.*t),2);
figure
plot(t, f(t), 'r', t, fFourier, 'b');
Warning: Imaginary parts of complex X and/or Y arguments ignored.
xlabel('t')
legend('f', 'Fourier approximation')
Bruno Luong
Bruno Luong 2023년 9월 9일
@Paul thanks
jvfa
jvfa 2023년 9월 10일

hi bruno, after studying your comment and code i finally understand. thank you, alsp for paul for the correction

댓글을 달려면 로그인하십시오.

추가 답변 (0개)

카테고리

도움말 센터File Exchange에서 Programming에 대해 자세히 알아보기

제품

릴리스

R2023a

질문:

2023년 9월 9일

편집:

2023년 9월 14일

Community Treasure Hunt

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

Start Hunting!

Translated by