terça-feira, 21 de abril de 2020

Como programar uma DFTS no Octave (ou Matlab)

Como programar uma DFTS no Octave (ou Matlab)

Das quatro representações de Fourier, a DTFS é a que é possível programar sem perdas nem na ida nem na volta pois não lidará com somas infinitas nem integrais, motivo pelo qual se torna interessante para se usar programando um computador para fazer os cálculos e é sobre isso que falaremos nesta postagem.

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.
    Como programar uma DFTS no Octave (ou Matlab)


    No Octave basta clicar no monitor preto entre o botão &gt; e a seta para esquerda azul.

    Como programar uma DFTS no Octave (ou Matlab)


    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.

    Como programar uma DFTS no Octave (ou Matlab)


    Como programar uma DFTS no Octave (ou 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

    Compartilhe:


    Algum link quebrado? Por favor, entre em contato para reportar esse erro.
    Leia a política de comentários do blog. Para escrever em $\LaTeX$ nos comentários, saiba em latex.mipedes.com.br.

    Postar um comentário

    Whatsapp Button works on Mobile Device only

    Digite no campo abaixo e tecle ENTER para pesquisar.