ft.mws

Fourierova transformace a její souvislost Fourierovou řadou
© Jiří Hospodka
Následující ukázka demonstroje použití programu Maple pro podporu výuky Fourierovy transformace v předmětech Elektrické obvody 2 (EO2) a počítačové řešení elektrických obvodů (PRO).

[Maple Bitmap]       [Maple Bitmap]  

Vyjděme z Fourierovy řady uvedeného periodického signálu a pokusme se odtud přejít k Fourierovu obrazu uvedeného obdélníkového pulzu. Většina operací je zatím bez komentáře (ty lze slyšet na cvičení :-).

Fourierova řada obdélníkového signálu podle obrázku, závistost spektrélních čar na parametech signálu  

>    f:=(t,t0::Range(0,0.5))->signum(cos(2*Pi*t)-cos(2*Pi*t0))/2+1/2;

f := proc (t, t0::Range(0,.5)) options operator, arrow; 1/2*signum(cos(2*Pi*t)-cos(2*Pi*t0))+1/2 end proc

>    rect:=plot(f(t,1/16),t=-0.5..2.5,color=blue):

>    rect;

[Maple Plot]

>    ck:=1/T*int(A*exp(-I*k*omega0*t),t=-tau/2..tau/2);

ck := -I/T*A*(exp(1/2*I*k*omega0*tau)-exp(-1/2*I*k*omega0*tau))/k/omega0

>    ck:=simplify(ck);

ck := 2/T*A*sin(1/2*k*omega0*tau)/k/omega0

>    ckn:=subs(A=1,tau=T/8,omega0=2*Pi/T,ck);

>    ck0:=limit(ckn,k=0);

ckn := sin(1/8*k*Pi)/k/Pi

ck0 := 1/8

>    sinx:=plot(1/8*sin(Pi*k/8)/(Pi*k/8),k=-20..20,color=blue):

>    sinx8:=plot([seq([[k,0],[k,ckn]],k=-20..-1),[[0,0],[0,ck0]],seq([[k,0],[k,ckn]],k=1..20)],color=red,thickness=3,xtickmarks=20):

>    plots[display](sinx,sinx8);

[Maple Plot]

>    limit(sin(x)/x,x=0);

1

Vztah ck  lze upravit na tvar

[Maple Bitmap]

  

x [Maple Bitmap]     

Kontrola

>    temp:=simplify(sum(ckn*exp(I*k*2*Pi*t),k=-10..-1)+ck0+sum(ckn*exp(I*k*2*Pi*t),k=1..10)):

>    pt:=plot(temp,t=-0.3..2.3,thickness=3,numpoints=200):

>    plots[display](pt,rect);

[Maple Plot]

>   

Změna periody T

Dosaďme nyní za tau = 1/8 , A = 1  a měňme periodu T od 1 do 5 s.

>    ckn:=subs(A=1,tau=1/8,omega0=2*Pi/T,ck):

>    ck0:=limit(ckn,k=0):

>    plots[animate](plot,[[seq([[k/T,0],[k/T,ckn]],k=-100..-1),[[0,0],[0,ck0]],seq([[k/T,0],[k/T,ckn]],k=1..100)],color=red,thickness=3,xtickmarks=20,view=[-20..20, 0.13..-0.03]],T=[1,2,3,4,5]);

[Maple Plot]

>   

Změna periody T pro součin ck*T

>    ck*T;

>    ckn:=subs(A=1,tau=1/8,omega0=2*Pi/T,T*ck):

>    ck0:=limit(ckn,k=0):

>    plots[animate](plot,[[seq([[k/T,0],[k/T,ckn]],k=-100..-1),[[0,0],[0,ck0]],seq([[k/T,0],[k/T,ckn]],k=1..100)],color=red,thickness=3,xtickmarks=20,view=[-20..20, 0.13..-0.03]],T=[1,2,3,4,5]);

2*A*sin(1/2*k*omega0*tau)/k/omega0

[Maple Plot]

>   

Změna šířky pulzu tau  pro součin ck*T   

>    cki:=subs(T=5,subs(A=1,tau=1/i,omega0=2*Pi/T,T*ck));

>    ck0:=limit(cki,k=0):

>    plots[animate](plot,[[seq([[k/5,0],[k/5,cki]],k=-100..-1),[[0,0],[0,ck0]],seq([[k/5,0],[k/5,cki]],k=1..100)],color=red,thickness=3,xtickmarks=20],i=[6,8,10,20,50,100]);

cki := 5*sin(1/5*k*Pi/i)/k/Pi

[Maple Plot]

>   

Změna šířky pulzu tau  pro součin ck*T při konstantní ploše pulzu tau *A

>    cki:=subs(T=5,subs(A=1*i,tau=1/i,omega0=2*Pi/T,T*ck));

>    ck0:=limit(cki,k=0):

>    plots[animate](plot,[[seq([[k/5,0],[k/5,cki]],k=-100..-1),[[0,0],[0,ck0]],seq([[k/5,0],[k/5,cki]],k=1..100)],color=red,thickness=3,xtickmarks=20],i=[6,8,10,20,50,100]);

cki := 5*i*sin(1/5*k*Pi/i)/k/Pi

[Maple Plot]

>   

Porovnejme výrazy pro koeficienty komplexní Fourierovy řady a Fourierovu transformaci

>    c[k] = ('1/T')*Int(f(t)*exp(-j*k*omega[0]*t),t = -T/2 .. T/2)

>    F(j*omega) = Int(f(t)*exp(-j*omega*t),t = -infinity .. infinity)

Zjistíme, že, ....

Exponenciální tvar

>    f(t) = Sum(c[k]*exp(j*k*omega[0]*t),k = -infinity .. infinity)

Zpětná Fourierova transformace

>    f(t) = ('1/(2*Pi)')*Int(F(j*omega)*exp(j*omega*t),omega = -infinity .. infinity)

Z předchozího je také evidentní, že Fourierův obraz Diracova pulzu je roven jeho ploše.

>    int(Dirac(t)*exp(-j*omega*t),t =-infinity..infinity);

Plocha Diracova pulzu.

>    int(Dirac(t),t=-infinity..infinity);

1

1

>   

Fourierova transformace pulzu

>    rect:=A*(Heaviside(t+tau/2)-Heaviside(t-tau/2));

rect := A*(Heaviside(t+1/2*tau)-Heaviside(t-1/2*tau))

>   

Jednotková funkce

>    plot(Heaviside(t),t=-1..1,thickness=2);

[Maple Plot]

>    plot(subs(A=1,tau=1/4,rect),t=-1..1,thickness=2);

[Maple Plot]

>    int(A*exp(-I*omega*t),t=-tau/2..tau/2);

-I*A*(exp(1/2*I*tau*omega)-exp(-1/2*I*tau*omega))/omega

>    F_rec:=simplify(%);

F_rec := 2*A/omega*sin(1/2*tau*omega)

>    with(inttrans);

[addtable, fourier, fouriercos, fouriersin, hankel, hilbert, invfourier, invhilbert, invlaplace, invmellin, laplace, mellin, savetable]

>    fourier(rect,t,omega);

2*A/omega*sin(1/2*tau*omega)

Jak bylo uvedeno, Fourierův obraz Diracova pulzu je roven jeho ploše. Nyní použijme pro výpočet přímý příkaz.

>    fourier(Dirac(t),t,omega);

1

>   

Nepoužívat funkci piecewise

>    rc:=piecewise(t>-tau/2 and t<tau/2,A);

rc := PIECEWISE([A, -1/2*tau < t and t < 1/2*tau],[0, otherwise])

>    plot(subs(A=1,tau=1/4,rc),t=-1..1,thickness=2);

[Maple Plot]

>    fourier(rc,t,omega);

fourier(PIECEWISE([A, -1/2*tau < t and t < 1/2*tau],[0, otherwise]),t,omega)

>   

Fourierova transformace periodických funkcí a vztah mezi Fourierovou tranzformací a řadou

Výsledkem Fourierova obrazu je tedy spojitá funkce, narozdíl od Fourierovy řady. Další odlišností je fyzikální rozměr této funkce (obrazu) oproti fyzikálnímu rozměru koeficientů řady (viz komentář dále). Pro náš případ obdélníhového puzlu dostaneme tedy velmi podobný výraz jeho obrazu v porovnání se vztahem pro fourierovy koeficienty periodicky se opakujícího průběhu. Stejný vztah dostaneme, pokud ve Fourierově obrazu přejdeme od spojitého kmitočtového spektra k diskrétnímu ( omega  -> k*omega0 ) a ten pak dělíme základní periodou signálu.

>    F_rec;

>    subs(omega=k*omega0,F_rec)/T;

>    ck;

2*A/omega*sin(1/2*tau*omega)

2/T*A*sin(1/2*k*omega0*tau)/k/omega0

2/T*A*sin(1/2*k*omega0*tau)/k/omega0

Obecně tedy platí následující pravidlo: Pokud známe Fourierův obraz časově omezené funkce s dobou trvání tau , můžeme jednoduše určit koeficienty Fourierovy řady periodické funkce, která vznikla z původní funkce jejím periodickým opakováním s periodou tau < T . Tyto koeficienty c[k]  určíme tak, že do Fourierova obrazu dosadíme omega = k*omega[0]  (přejdeme od spojitého kmitočtového spektra k diskrétnímu, kde omega[0] = 2*Pi/T ) a výsledek pak dělíme periodou T . To také vyplývá z výše uvedených vztahů - ze vztahu pro koeficienty Fourierovy řady c[k]  a pro Fourierův obraz F(j*omega) .

Použijme toto pravidlo na náš případ. Mějme symetrický pulz o šířce tau = 1/2  . ...

>    plot(subs(A=1,tau=1/2,rect),t=-1..1,thickness=2);

>    F_obd:=fourier(rect,t,omega);

[Maple Plot]

F_obd := 2*A/omega*sin(1/2*tau*omega)

>    F_obd:=fourier(rect,t,omega);

F_obd := 2*A*sin(1/2*omega*tau)/omega

>    plot(signum(cos(2*Pi*t)),t=0..1.5,thickness=2);

>    a_obd:=8/T*int(1*cos(k*2*Pi/T*t), t=0..T/4 );

[Maple Plot]

a_obd := 4*sin(1/2*k*Pi)/k/Pi

>    ck:=subs(A=2,omega=k*2*Pi/(tau*2),F_rec)/(2*tau);

ck := 2*sin(1/2*k*Pi)/k/Pi

>    c0:=limit(ck,k=0);

c0 := 1

>    evalc(c0+sum(ck*exp(I*2*Pi*k*t),k=-10..-1)+sum(ck*exp(I*2*Pi*k*t),k=1..10));

1+4/9*1/Pi*cos(18*Pi*t)-4/7*1/Pi*cos(14*Pi*t)+4/5*1/Pi*cos(10*Pi*t)-4/3*1/Pi*cos(6*Pi*t)+4/Pi*cos(2*Pi*t)

>    sum(a_obd*cos(2*Pi*k*t),k=1..10);

4/Pi*cos(2*Pi*t)-4/3*1/Pi*cos(6*Pi*t)+4/5*1/Pi*cos(10*Pi*t)-4/7*1/Pi*cos(14*Pi*t)+4/9*1/Pi*cos(18*Pi*t)

>   

>    fourier(A*exp(I*omega0*t),t,omega) assuming omega0>0;

>    fourier(A*cos(omega0*t),t,omega) assuming omega0>0;

>    fourier(A*sin(omega0*t),t,omega) assuming omega0>0;

2*A*Pi*Dirac(omega-omega0)

A*Pi*(Dirac(-omega+omega0)+Dirac(omega+omega0))

A*Pi*(-Dirac(omega-omega0)+Dirac(omega+omega0))*I

Tyto vztahy lze poměrně snadno odvodit při znalosti spektra konstantního signálu ( [Maple Bitmap] , což lzde dokázat zpětnou fransformací [Maple Bitmap]  ) a použitím věty o kmitočtovém posunu

[Maple Bitmap]

 Lze tak získat obraz komplexní exponencialy.

[Maple Bitmap]   což lze potvrdit zpětnou transformací   [Maple Bitmap]  

Odtud již lze napsat i ostatní vztahy (pro funkce sin a cos), pokud uvážíme, že [Maple Bitmap] . Podobně odvodit i vztah pro obecný periodický signál

[Maple Bitmap]   kde   [Maple Bitmap]  

Jaký je tedy vztah mezi spektrem periodického signálu vypočteným pomocí Fourierovy transformace a koeficienty jeho Fourierova rozvoje. Například pro spektrum funkce A*cos(omega0*t)  platí výše uvedený vztah. Vzhledem k~tomu, že A*cos(omega0*t) = A/2*(exp(I*omega0*t)+exp(-I*omega0*t)) , jsou koeficienty komplexní Fourierovy řady c[-1] = c[1]  = A/2 . Jejich grafické znázornění je uvedeno na následijícím obrázku (a). Oproti tomu spektrum stejného signálu A*cos(omega0*t) , získané Fourierovou transformací tvoří dva Diracovy pulzy s plochou Pi*A  na kmitočtech -omega  a omega , jak ukazuje obrázek (b).

[Maple Bitmap]

Plocha Diracových pulzů je přitom úměrná velikosti koeficientu ck  Fourierova rozvoje. Je-li fyzikální rozměr amplitudy A  např. [V], pak plocha tohoto pulzu má fyzikální rozměr [V * s * rad / s]=[V * rad]. Pokud tedy chceme znát velikost příslušného koeficientu ck  stačí vydělit plochu odpovídajícího Diracova pulzu členem 2*Pi . To také evidentní ze vztahu pro zpětnou Fourierovu transformaci. Pokud ho budeme aplikovat na dané spektrum (např. kosinusového signálu), vynásobené amplitudou A , dostaneme člen 1/(2*Pi)  násobený integrálem, udávajícím plochy obou pulzů násobené komplexními exponenciálami, tj.

[Maple Bitmap] [Maple Bitmap]

Po úpravě samozřejmě dostaneme původní signál se správnou velikostí amplitudy.

>   

Fourierova transformace trojůhelníkového pulzu

>    troj:=piecewise(t>-t0 and t<=0,(t+t0)/t0*U,t>0 and t<t0,(t0-t)/t0*U,0);

troj := PIECEWISE([(t+t0)/t0*U, -t0 < t and t <= 0],[(t0-t)/t0*U, 0 < t and t < t0],[0, otherwise])

>    plot(subs(t0=1,U=2,troj),t=-2..2,thickness=2);

[Maple Plot]

>    F:=int(troj*exp(-I*omega*t),t=-infinity..infinity) assuming t0>0;

F := U*(-exp(omega*t0*I)+1+omega*t0*I)/omega^2/t0-U*(exp(omega*t0*I)*omega*t0*I-exp(omega*t0*I)+1)*exp(-I*t0*omega)/omega^2/t0

>    Fu:=simplify(F);

Fu := -2*U*(cos(omega*t0)-1)/omega^2/t0

>    fourier(subs(t0=1,U=2,troj),t,omega);combine(%);

8/omega^2*sin(1/2*omega)^2

(4-4*cos(omega))/omega^2

>    simplify(subs(t0=1,U=2,Fu));

-4*(-1+cos(omega))/omega^2

>    plot(subs(t0=1,U=2,Fu),omega=-4*Pi..4*Pi,thickness=2);

[Maple Plot]

>    ft:=invfourier(Fu,omega,t) assuming t0>0;

ft := U/t0*((t-t0)*Heaviside(t-t0)+(t+t0)*Heaviside(t+t0)-2*t*Heaviside(t))

>    convert(ft,piecewise,t) assuming t0>0;

PIECEWISE([0, t < -t0],[U/t0*undefined, t = -t0],[(t+t0)/t0*U, t < 0],[U/t0*(t0+undefined), t = 0],[(t0-t)/t0*U, t < t0],[U/t0*undefined, t = t0],[0, t0 < t])

>    laplace(subs(t0=1,U=2,troj),t,p);

2/p-2/p^2*(1-exp(-p))

>   

>   

>   

Laplaceova transformace

>    F(p) = Int(f(t)*exp(-p*t),t = 0 .. infinity)

Zpětná Laplaceova transformace

>    f(t) = ('1/(2*Pi*j)')*Int(F(p)*exp(p*t),p = sigma-j*infinity .. sigma+j*infinity)

>   

Úloha 2.1.1

>    ex:=U*exp(-a*t);

ex := U*exp(-a*t)

>    plot(subs(U=1,a=3,ex)*Heaviside(t),t=-1..2,thickness=2);

[Maple Plot]

>    int(ex*exp(-I*omega*t),t=-infinity..infinity) assuming a>0;

undefined

>    int(ex*exp(-I*omega*t),t=0..infinity) assuming a>0;

1/(a+omega*I)*U

>    int(ex*Heaviside(t)*exp(-I*omega*t),t=-infinity..infinity) assuming a>0;

1/(a+omega*I)*U

Pozor, jednotkový skok je nutný.

>    fourier(ex,t,omega) assuming a>0;

>    fourier(ex*Heaviside(t),t,omega) assuming a>0;

U*fourier(exp(-a*t),t,omega)

1/(a+omega*I)*U

>    ex2:=U*exp(-a*(t-t0));

ex2 := U*exp(-a*(t-t0))

>    plot(subs(U=1,a=3,t0=0.2,ex2*Heaviside(t-t0)),t=-1..2,thickness=2);

[Maple Plot]

>    fourier(ex2*Heaviside(t-t0),t,omega) assuming a>0, t0>0;

U/(a+omega*I)*exp(-I*omega*t0)

>   

Úloha 2.1.3

Máme určit frekvenční přenos obvodu (do kmitočtu fm , kde už je přenos obvod u téměř nulový. Obvod máme vybudit pulzem a kmitočtový přenos pak určit z jeho odezvy. Jak úzký pulz je třeba pro vybuzení obvodu?

Spektrum pulzu.

>    F_rec;

2*A*sin(1/2*omega*tau)/omega

Lépe je to vidět z tohoto vztahu [Maple Bitmap] . Dosaďme omega = 2*Pi*f , tj. omega*tau/2 = Pi*f*tau  a označme součin f*t = x . Pak vyneseme, jak se zmenšuje spektrum v závislosti na tomto součinu (relativně vůči ploše A*tau ).  

>    plot(sin(Pi*x)/(Pi*x),x=0..1);

[Maple Plot]

Jestliže zvolíme fm*tau = 1/10 , pak spektrum takového pulzu vykazuje pro tento kmitočet f = fm  pokles cca o 2%, což je možné tolerovat.

>    evalf(subs(x=0.1,sin(Pi*x)/(Pi*x)));

.9836316430

Prakticky např. pro fm = 100 kHz vychází tau = 1*mu s.

>    solve(100e3*tau=0.1);

.1000000000e-5

>   

Úloha 2.1.4

>    S:=1/1000;

>    u2:=1/2*exp(-1000*t)*Heaviside(t);

S := 1/1000

u2 := 1/2*exp(-1000*t)*Heaviside(t)

>    U1:=fourier(Dirac(t)*S,t,omega);

>    U2:=fourier(u2,t,omega);

U1 := 1/1000

U2 := 1/(2*(1000+omega*I))

>    P:=U2/U1;

P := 500/(1000+omega*I)

>    evalf(numer(P)/1000,1)/evalf((denom(P)/1000),1);

.5/(1.+.1e-2*I*omega)

>   

Příklad 2.1.4

>    u1:=2*exp(-2000*t)*Heaviside(t);

>    u2:=(8*exp(-2000*t)-6*exp(-3000*t)-2*exp(-1000*t))*Heaviside(t);

u1 := 2*exp(-2000*t)*Heaviside(t)

u2 := (8*exp(-2000*t)-6*exp(-3000*t)-2*exp(-1000*t))*Heaviside(t)

>    U1:=fourier(u1,t,omega);

>    U2:=fourier(u2,t,omega);

>    P:=U2/U1;

U1 := 2/(2000+omega*I)

U2 := 4000*I*omega/(2000+omega*I)/(3000+omega*I)/(1000+omega*I)

P := 2000*I*omega/(3000+omega*I)/(1000+omega*I)

>   

Výpočty pomocí Laplaceovy transformace

Laplaceova transformace Diracova pulzu.

>    laplace(Dirac(t),t,p);

1

Laplaceova transformace jednotkového skoku.

>    laplace(Heaviside(t),t,p);

1/p

Příkaz "laplace" evidentně ignoruje funkční hodnoty pro t < 0  (integruje se od 0).

>    laplace(1,t,p);

1/p

Následují výpočty některých výše uvedených příkladů pomocí Laplaceovy transformace.

>    laplace(U*exp(-a*t),t,p);

U/(p+a)

Je nutné vynásobit jednotkovým skokem.

>    laplace(U*exp(-a*(t-t0)),t,p);

>    laplace(U*exp(-a*(t-t0))*Heaviside(t-t0),t,p) assuming t0>0;

U*exp(a*t0)/(p+a)

U/(p+a)*exp(-p*t0)

>    invlaplace(%,p,t) assuming  t0>0;

U*exp(-a*(t-t0))*Heaviside(t-t0)

>    u1:=2*exp(-2000*t)*Heaviside(t);

>    u2:=(8*exp(-2000*t)-6*exp(-3000*t)-2*exp(-1000*t))*Heaviside(t);

u1 := 2*exp(-2000*t)*Heaviside(t)

u2 := (8*exp(-2000*t)-6*exp(-3000*t)-2*exp(-1000*t))*Heaviside(t)

>    U1:=laplace(u1,t,p);

>    U2:=laplace(u2,t,p);

>    P:=U2/U1;

U1 := 2/(p+2000)

U2 := 4000*p/(p+2000)/(p+3000)/(p+1000)

P := 2000*p/(p+3000)/(p+1000)

Předpokládá se, že funkce "začíná" v nule.

>    invlaplace(U2,p,t);

8*exp(-2000*t)-6*exp(-3000*t)-2*exp(-1000*t)

>