function [xs, g, PSDx, PSDy] = wiener_fd(x, y, sigma_noise)
    N = length(x);
    Y = fft(y);
    X = fft(x);
    PSDx = 1/N*abs(X).^2; % spectral power
    PSDy = 1/N*abs(Y).^2; % spectral power

    k = sqrt(sum(PSDx)./sum(PSDy))
    % fourier transform of the wiener filter
    G = PSDx./(PSDx + sigma_noise.^2); 
    % perform the filtering over the fourier
    G = 1/k*G.*PSDx./PSDy;
    xs = ifft(G.*Y);
    g = ifft(G);