Transformada discreta de Fourier y su inversa usando MATLAB

Con la llegada de MATLAB y toda la incorporación científica que ha traído, ha habido un cambio significativo y una simplificación de escenarios de ingeniería sofisticados. De hecho, este cambio ha contribuido a ayudar a desarrollar una mejor visualización y habilidades prácticas para los estudiantes que cursan estudios técnicos en ciencias, si no en otros campos al menos. Aquí analizamos la implementación de una idea matemática fundamental: la transformada discreta de Fourier y su inversa utilizando MATLAB. 

Cálculo de la DFT

Las ecuaciones estándar que definen cómo la transformada discreta de Fourier y la inversa convierten una señal del dominio del tiempo al dominio de la frecuencia y viceversa son las siguientes: 

DFT:  para k=0, 1, 2….., N-1

IDFT:  para n=0, 1, 2….., N-1

Siendo las ecuaciones bastante sencillas, uno podría simplemente ejecutar bucles repetitivos/anidados para la suma y terminar con eso. Sin embargo, debemos intentar utilizar otro método en el que usemos arrays para encontrar la solución al problema. Muchos lectores recordarán que la DFT y la IDFT de una señal en el dominio del tiempo/frecuencia se pueden representar en formato vectorial de la siguiente manera:

Cuando tomamos los factores de twiddle como componentes de una array, se vuelve mucho más fácil calcular la DFT y la IDFT. Por lo tanto, si nuestra señal en el dominio de la frecuencia es una array de una sola fila representada por X _N y la señal en el dominio del tiempo también es una array de una sola fila representada como x _N ……

Con esta interpretación, todo lo que necesitamos hacer es crear dos arrays sobre las cuales emitiremos una multiplicación de arrays para obtener la salida. La array de salida SIEMPRE será una array de orden Nx1 ya que tomamos una array de una sola fila como nuestra señal de entrada (X _N o x _N ). Este es esencialmente un vector que podemos transponer a una array horizontal para nuestra conveniencia.

Algoritmo (DFT)

1. Obtenga la secuencia de entrada y el número de puntos de la secuencia DFT.
2. Envía los datos obtenidos a una función que calcula el DFT. No es imperativo declarar una nueva función, pero la legibilidad y el flujo del código se vuelven más claros y evidentes.
3. Determine la longitud de la secuencia de entrada utilizando la función length( ) y verifique si la longitud es mayor que el número de puntos. N siempre debe ser igual o mayor que la secuencia. Si intenta ejecutar la multiplicación de arrays sin cumplir la condición, aparecerá un error en la ventana de comandos.
4. Contabilización de la diferencia en longitudes de la secuencia de entrada y los N puntos usando una array separada que agrega ceros adicionales para alargar la secuencia de entrada. Esto se hace usando la función zeros(no_of_rows, no_of_columns) que crea una array 2D compuesta de ceros.
5. Con base en el valor de N obtenido como entrada, cree la array W _N. Para hacer esto, implemente 2 bucles ‘for’, un procedimiento bastante básico.
6. Simplemente multiplique las dos arrays que se han creado. Esta es una array de las muestras de señal de dominio de frecuencia requeridas.
7. Trace la magnitud y la fase de la señal de salida a través de las funciones incorporadas abs(function_name) y angle(function_name) .

% MATLAB code for DFT 
clc;
xn=input('Input sequence: ');
N = input('Enter the number of points: ');
Xk=calcdft(xn,N);
disp('DFT X(k): ');
disp(Xk);
mgXk = abs(Xk);
phaseXk = angle(Xk);
k=0:N-1;
subplot (2,1,1);
stem(k,mgXk);
title ('DFT sequence: ');
xlabel('Frequency');
ylabel('Magnitude');
subplot(2,1,2);
stem(k,phaseXk);
title('Phase of the DFT sequence');
xlabel('Frequency');
ylabel('Phase');
  
function[Xk] = calcdft(xn,N)
    L=length(xn);
    if(N<L)
        error('N must be greater than or equal to L!!')
    end
    x1=[xn, zeros(1,N-L)];
    for k=0:1:N-1
        for n=0:1:N-1
            p=exp(-i*2*pi*n*k/N);
            W(k+1,n+1)=p;
        end
    end
    disp('Transformation matrix for DFT')
    disp(W);
    Xk=W*(x1.')
end

Producción:

>> Input sequence: [1 4 9 16 25 36 49 64 81]
>> Enter the number of points: 9

Algoritmo (IDFT)

1. Obtenga la señal/secuencia en el dominio de la frecuencia como entrada ( X(k) ). La longitud de esta secuencia es suficiente como valor para N (puntos).
2. Pase esta array a una función para el cálculo.
3. Ejecute 2 bucles en la función para crear la array. Tenga en cuenta que esta array debe conjugarse cuando se utiliza para el cálculo. Puede optar por declarar explícitamente otra array que sea el conjugado de la array W _N .
4. Una vez que se ha creado la array, obtenga el conjugado usando ‘ * ‘ y simplemente multiplíquelo con la transposición de la secuencia de entrada. Requerimos la transposición ya que la entrada es una array de fila. Al multiplicar con la array W _N que hemos creado, el número de columnas en W _N debe coincidir con el número de filas en X(k).
5. Traza esta secuencia usando stem(x_axis, y_axis) . NO use plot( ) ya que esta no es una señal CT.

% MATLAB code for IDFT
clc;
Xk = input('Input sequence X(k): ');
xn=calcidft(Xk);
N=length(xn);
disp('xn');
disp(xn);
n=0:N-1;
stem(n,xn);
xlabel('time');
ylabel('Amplitude');
  
function [xn] = calcidft(Xk) %function to calculate IDFT
    N=length(Xk);
    for k=0:1:N-1
        for n=0:1:N-1
            p=exp(i*2*pi*n*k/N);
            IT(k+1,n+1)=p;
        end
    end
    disp('Transformation Matrix for IDFT');
    disp(IT);
    xn = (IT*(Xk.'))/N;
end

Producción:

>> Enter the input sequence: [1 2 3 4 5 9 8 7 6 5]

Array de transformación

La secuencia en el dominio del tiempo