Referencial Teórico
Considere um sinal x[n] com período fundamental N (natural) então temos$$x[n]=\sum_{k\in \left\langle N \right\rangle} X[k]e^{jk\Omega_0 n}$$ $$X[k]=\frac{1}{N}\sum_{n\in \left\langle N \right\rangle} x[n]e^{-jk\Omega_0 n}$$ onde $\Omega_0=\frac{2\pi}{N}$ e $X[k]$ são os coeficientes de Fourier (da forma complexa), $j^2=-1$ e $$e^{j\Box}=\cos(\Box)+j.\sin(\Box)$$.
A relação entre $x[n]$ e $X[k]$ permite a representação de um sinal em domínios diferentes. O primeiro no que pode ser chamado de domínio do tempo e o segundo no domínio da frequência. Escrever um programa de computador que faça esse trabalho é interessante do ponto de vista pedagógico pois coloca o estudante para pensar no que está sendo executado, tanto para a ida quanto para a volta. Aqui falaremos da ida, ou seja, sair de um sinal no domínio do tempo e ir para o domínio da frequência mostrando uma saída gráfica de $|X[k]|$ e $Arg\{X[k]\}$, os espectros de magnitude e fase, respectivamente. Em outro momento podemos escrever o programa da volta, ou seja, trazer para o domínio do tempo um sinal representado no domínio da frequência por seus espectros de magnitude e fase.
Uma observação rápida: DTFS não é a forma mais rápida de fazer essa transformação. FFT é mais rápida, mas o objetivo aqui é programar o computador através do Octave ou Matlab para executar essa ação. Essas funções já existem nos respectivos programas. Não se trata de nada novo aqui, mas só de treino mesmo. ;-)
Como sua função será chamada?
O código mostrado abaixo pode ser colocado em um arquivo *.m que rodará tanto no Octave quanto no MATLAB. Ele mostrará os Espectros de Magnitude e de Fase de um sinal x[n] com período N. O procedimento mostra ao final uma janela com os referidos espectros no intervalo [-N+1,N], ou seja, colocarmos o sinal com dois períodos fundamentais para que ficasse mais fácil a visualização.Para chamar a função basta escrever DTFS( <vetor>). Por exemplo:
- n=1:21
- xn=cos(2*pi/21*n)
- DTFS(xn)
Saídas pelo Octave e pelo Matlab
O que o procedimento faz é simples. Você entra com o sinal x[n] e o tamanho do período ele mostra X[k]. Como X[k] não é, necessariamente, real, mostramos o |X[k]| (Espetro de Magnitude) e o Arg{X[k]} (Espectro de Fase).Como mencionado anteriormente é possível rodar o arquivo *.m tanto no Octave quanto no MATLAB.
No Octave basta clicar no monitor preto entre o botão > e a seta para esquerda azul.
e no MATLAB (figura anterior) basta apertar a tecla F5.
Em ambos a resposta será como mostrado nas figuras seguintes, respectivamente, para o Octave e o MATLAB.
O código que geral esta saída está em azul logo a seguir. Naturalmente que % serve para comentar e que, naturalmente, o código pode ser melhorado.
function DTFS(x) N=size(x,2); Omega0=2*pi/N; Xr=zeros(1,N); Xi=zeros(1,N); EspMag=zeros(1,N); EspFase=zeros(1,N); for k=1:N for n=1:N Xr(k)=Xr(k)+1/N*cos(k*Omega0*n)*x(n); Xi(k)=Xi(k)-1/N*sin(k*Omega0*n)*x(n); end end for k=1:N if abs(Xr(k))<0.0001 Xr(k)=0; end if abs(Xi(k))<0.0001 Xi(k)=0; end end for k=1:N if (Xr(k)>0) EspFase(k)=atan(Xi(k)/Xr(k)); elseif (Xr(k)<0 && Xi(k)>0) EspFase(k)=atan(Xi(k)/Xr(k))+pi; elseif (Xr(k)<0 && Xi(k)<0) EspFase(k)=atan(Xi(k)/Xr(k))-pi; elseif (Xr(k)==0)&&(Xi(k)>0) EspFase(k)=pi/2; elseif (Xr(k)==0)&&(Xi(k)<0) EspFase(k)=-pi/2; else EspFase(k)=0; end end for k=1:N EspMag(k)=sqrt(Xr(k)^2+Xi(k)^2); end Xr2=zeros(1,2*N); Xi2=zeros(1,2*N); for k=1:N Xr2(k)=EspMag(k); Xi2(k)=EspFase(k); end for k=N+1:2*N Xr2(k)=EspMag(k-N); Xi2(k)=EspFase(k-N); end subplot(2,1,1); stem(-N+1:N,Xr2) xlabel('k'); ylabel('|X(k)|'); subplot(2,1,2); stem(-N+1:N,Xi2) xlabel('k'); ylabel('Arg(X(k))');
Grande abraço
Luís Cláudio LA
Postar um comentário