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