[matlab] Discrete Fourier Transformation (DFT) (이산 퓨리에 변환!!!, 필터, 신호 및 시스템, 익명 함수, for, Signals and Systems)

0. [matlab] - Discrete Fourier Transformation (DFT)

1. 풀이

이번 신호 및 시스템을 달려오며 최종적인 목표인 DFT를 수행해보았다.

이산시간 퓨리에 변환은 행렬 곱셈으로 DFT를 수행하도록 한 것이다!!!

이런 방식을 통해 컴퓨터에서 빠르게 DFT를 수행하여 원하는 값을 얻을 수 있고, 또한 이 변환은 역변환 또한 행렬 곱을 통해 손쉽게 얻을 수 있다.

이렇게 얻은 퓨리에 변환을 활용하는 과제가 2번째 과제였는데, 나의 코드는 사실 올바른 필터의 구현이 아니다.

이미 DFT 변환이 완료된 형태이기 때문에, 필터와 convolution을 수행하는 것이 아닌, 행렬 전체에 3을 곱하면 원하는 결과을 얻는 것이다.

이렇게 수행하면, 이미지가 전체적으로 밝아지는 결과를 얻을 수 있을 것이다.

2. 소스코드

1.

clear;
 
%% (iii) 사진에 DFT 적용
I = imread('Lenna_New.png');            % 사진 데이터 행렬을 I에 옮긴다.
I_gray = rgb2gray(I);                   % 사진이 3차원이므로 2차원으로 변환한다.(무채색으로 변환)
I_gray_double = double(I_gray);         % 데이터를 double형으로 변환
 
X_N = my_DFT2(I_gray_double);           % two-dimensional DFT 수행
 
m = abs(X_N);                           % Megnitude
% X_N(m<1e-6) = 0;                      % 크기가 작은 변환 값들을 0으로 설정하여 반올림
                                        % 오차를 줄인다.(하지만 생략해주었다.)
                                            
p = unwrap(angle(X_N));                 % unwrap를 통해 Phase의 불연속점을 제거한다.
f = (0:length(X_N)-1)*100/length(X_N);  % Frequency vector
 
% 결과 출력
figure(1);
subplot(2,1,1)
plot(f,m)
title('Magnitude')
 
subplot(2,1,2)
plot(f,p*180/pi)
title('Phase')
 
 
%% (iv) fft2 활용 DFT
X_MATLAB = fft2(I_gray_double);
 
m = abs(X_MATLAB);                               % Megnitude
% y(m<1e-6) = 0;
p = unwrap(angle(X_MATLAB));                       % Phase
f = (0:length(X_MATLAB)-1)*100/length(X_MATLAB);      % Frequency vector
 
figure(2);
subplot(2,1,1)
plot(f,m)
title('Magnitude')
 
subplot(2,1,2)
plot(f,p*180/pi)
title('Phase')
 
 
%% (v) 오차 계산
E = X_N - X_MATLAB;
det(E'*E)


%% (vi) 스케일링된 색으로 이미지 표시
% 기본 적용
figure(3);
subplot(2,1,1); 
imagesc(abs(fftshift(X_N))); colorbar; title('my DFT');
subplot(2,1,2);
imagesc(abs(fftshift(X_MATLAB))); colorbar; title('fft2');

% log를 활용해 더 자세하게 분석
figure(4);
subplot(2,1,1); 
imagesc(abs(log2(fftshift(X_N)))); colorbar; title('my DFT log');
subplot(2,1,2);
imagesc(abs(log2(fftshift(X_MATLAB)))); colorbar; title('fft2 log');


%% (vii) 다시 그림으로 변환
x_IDFT = my_inv_DFT2(X_N);

% 다시 int형으로 반환하여 그림으로 다시 출력된 것을 test
 I_IDFT = uint8(x_IDFT)
 figure(5);
 imshow(I_IDFT);


%% (viii) 오차 계산
F = x_IDFT-I_gray_double;
det(F'*F);
 
 
%% (i) (ii) two-dimensional DFT
function X_N = my_DFT2(x)
    %% (i)번 행렬 구하기
    [x_row, x_col] = size(x);       % x 행렬의 크기
    W_M = zeros(x_row, x_row);      % W_M 행렬 [x row*x row]
    W_N = zeros(x_col, x_col);      % W_N 행렬 [x col*x col]
    
    % 행렬을 알고리즘에 맞게 채워준다.
    for a = 1 : 1 : x_row
        for b = 1 : 1 : x_row
            W_M(a,b) = exp((a-1)*(b-1)*-j*2*pi/x_row);
        end
    end
    
    for a = 1 : 1 : x_col
        for b = 1 : 1 : x_col
            W_N(a,b) = exp((a-1)*(b-1)*-j*2*pi/x_col);
        end
    end
    %%  (ii)번 DFT 수행
    % DFT 변환한 X를 반환한다.
    X_N = W_M*x*W_N;
end
 
%% (vii) inverse DFT 수행
function x = my_inv_DFT2(X)
    [x_row, x_col] = size(X);
    W_M = zeros(x_row, x_row);
    W_N = zeros(x_col, x_col);
 
    for a = 1 : 1 : x_row
        for b = 1 : 1 : x_row
            W_M(a,b) = exp((a-1)*(b-1)*-j*2*pi/x_row);
        end
    end
    
    for a = 1 : 1 : x_col
        for b = 1 : 1 : x_col
            W_N(a,b) = exp((a-1)*(b-1)*-j*2*pi/x_col);
        end
    end
    x = inv(W_M)*X*inv(W_N);
end
 
 

2.

clear;
 
%% (i) 사진에 DFT 적용
I = imread('Lenna_New.png');    % 사진 데이터 행렬을 I에 옮긴다.
I_gray = rgb2gray(I);           % 사진이 3차원이므로 2차원으로 변환한다.(무채색으로 변환)
I_gray_double = double(I_gray);         % 데이터를 double형으로 변환
 
X_N = my_DFT2(I_gray_double);   % two-dimensional DFT 수행
 
%% (ii) DFT 변환된 X행렬에 filter적용
H = ones(5,4)*3;
% X_filter_after = my_DFT2_filter(X_N, H);
X_filter_after = filter2(H,X_N,'same');
X_filter_after =3*X_N;
%% (iii) 다시 그림으로 변환
x_IDFT = my_inv_DFT2(X_filter_after);
 
% 그림 출력
I_IDFT = uint8(x_IDFT);
figure(1);
imshow(I_IDFT);
 
%%%%%%
% imwrite(I_IDFT,'Lenna_filter.png')        %파일 입출력
%%%%%%
 
%% 아래는 함수 구현 파트
 
%% two-dimensional DFT
function X_N = my_DFT2(x)
    [x_row, x_col] = size(x);       % x 행렬의 크기
    W_M = zeros(x_row, x_row);      % W_M 행렬 [x row*x row]
    W_N = zeros(x_col, x_col);      % W_N 행렬 [x col*x col]
    
    % 행렬을 알고리즘에 맞게 채워준다.
    for a = 1 : 1 : x_row
        for b = 1 : 1 : x_row
            W_M(a,b) = exp((a-1)*(b-1)*-j*2*pi/x_row);
        end
    end
    
    for a = 1 : 1 : x_col
        for b = 1 : 1 : x_col
            W_N(a,b) = exp((a-1)*(b-1)*-j*2*pi/x_col);
        end
    end
    % DFT 변환한 X를 반환한다.
    X_N = W_M*x*W_N;
end
 
%% inverse DFT
function x = my_inv_DFT2(X)
    [x_row, x_col] = size(X);
    W_M = zeros(x_row, x_row);
    W_N = zeros(x_col, x_col);
 
    for a = 1 : 1 : x_row
        for b = 1 : 1 : x_row
            W_M(a,b) = exp((a-1)*(b-1)*-j*2*pi/x_row);
        end
    end
    
    for a = 1 : 1 : x_col
        for b = 1 : 1 : x_col
            W_N(a,b) = exp((a-1)*(b-1)*-j*2*pi/x_col);
        end
    end
    x = inv(W_M)*X*inv(W_N);
end
 
%% 2-dimension convolution을 활용해 filter를 구현했다.
function con = my_DFT2_filter(x, h)
    [x_row, x_col] = size(x);
    [h_row, h_col] = size(h);
    con = zeros(x_row, x_col);
    
    % 이중 for문을 활용해 일정 영역 내부의 x*h를 전부 더한다.
    for m = 1:1:x_row
        for n = 1:1:x_col
            for u=1:1:h_row
                for v=1:1:h_col
                    if (m+u<x_row+1) && (n+v<x_col+1)
                        con(m,n) = con(m,n) + h(u,v)*x(m+u-1,n+v-1);
                    end
                end
            end
        end
    end
end
 
 
 

3. 참고

질문이나 지적 있으시면 댓글로 남겨주세요~

도움 되셨으면 하트 꾹!