nf2ff.m

a usefull source code from internet

% Near field to far field transformation for near field 
% antenna measurement setups
% K. Van Caekenberghe, M. J. Pelk
% 
% The code uses:
% 1. GoldsteinUnwrap2D.m, posted on the webstie by B. Spottiswoode on 
%    22 December 2008
% 2. Sphere3D.m, posted on the fileexchange website by J. M. De Freitas on 
%    15 September 2005
% 
% Sought-after extensions are listed below. Please post them on the 
% website, if you are willing to contribute:
% 1. Cross polarization (Ludwig 1, 2, and 3)
% 2. Probe compensation
% 3. Time gating
%
% Please report bugs by e-mail ([e-mail], [e-mail]) or
% through the website.
%
% References:
% 1. C. A. Balanis, "Antenna Theory, Analysis and Design, 2nd Ed.", Wiley, 
%    1997.
% 2. R. M. Goldstein, H. A. Zebken, and C. L. Werner, "Satellite Radar 
%    Interferometry: Two-Dimensional Phase Unwrapping", Radio Sci., 
%    vol. 23, no. 4, pp. 713-720, 1988.
% 3. D. C. Ghiglia and M. D. Pritt, "Two-Dimensional Phase Unwrapping:
%    Theory, Algorithms and Software". Wiley-Interscience, 1998.
% 4. J. M. De Freitas. "SPHERE3D: A Matlab Function to Plot 3-Dimensional 
%    Data on a Spherical Surface". 
%    QinetiQ Ltd, Winfrith Technology Centre, Winfrith, 
%    Dorchester DT2 8XJ. UK. 15 September 2005. 

clc;
close all;
clear all;

c=299792458; % Speed of light in vacuum [m/s]
load('scanarray_pol2_h6mm-10_2_2009.mat'); % Measurement data
sdata2=sdata;
load('scanarray_pol1_h6mm-10_2_2009.mat'); % Measurement data
freq=sdata.freq; % Measured frequency points [Hz]
N=length(freq);
f_start=freq(1); 
f_stop=freq(N);
df=(f_stop-f_start)/(N-1);
dt=1/(N*df);
t=(0:N-1)*dt;
x=c*t;

figure;
subplot(2,1,1);
plot(freq/1e9,20*log10(abs(sdata.s21{floor(sdata.ypoints/2),floor(sdata.xpoints/2)})),'k');
xlabel('Frequency (GHz)');
ylabel('S_{21} (dB)');
subplot(2,1,2);
plot(freq/1e9,180/pi*angle(sdata.s21{floor(sdata.ypoints/2),floor(sdata.xpoints/2)}),'k');
xlabel('Frequency (GHz)');
ylabel('\angle S_{21} (°)');

figure;
subplot(2,1,1);
plot(t'*1e9,20*log10(abs(ifft(squeeze(sdata.s21{floor(sdata.ypoints/2),floor(sdata.xpoints/2)})))),'k');
xlabel('Time (ns)');
ylabel('S_{21} (dB)');
subplot(2,1,2);
plot(x',20*log10(abs(ifft(squeeze(sdata.s21{floor(sdata.ypoints/2),floor(sdata.xpoints/2)})))),'k');
xlabel('Distance (m)');
ylabel('S_{21} (dB)');

clear N, x;

theta=[-pi/2+0.1:0.05:pi/2-0.1];
phi=[0+0.1:0.05:pi-0.1];
[theta, phi]=meshgrid(theta,phi);

% See equations (16-10a) and (16-10b) in Balanis
M = sdata.xpoints; % Amount of samples in the x direction (along table, left to right)
N = sdata.ypoints; % Amount of samples in the y direction (across table, front to back)
dx = sdata.x_step/1000; % Sample spacing in the x direction [m]
dy = sdata.y_step/1000; % Sample spacing in the y direction [m]
a = dx*(M-1); % The length of the scanned area in the x direction [m]
b = dy*(N-1); % The length of the scanned area in the y direction [m]
x = [-a/2:a/(M-1):a/2];
y = [-b/2:b/(N-1):b/2];

% Zero padding is used to increase the resolution of the plane wave spectral domain.
MI=2^(ceil(log2(M))+1);
NI=2^(ceil(log2(N))+1);

% See equations (16-13a) and (16-13b) in Balanis
% Substitute M and N with MI and NI respectively
m=[-MI/2:1:MI/2-1];
n=[-NI/2:1:NI/2-1];

k_X_Rectangular=2*pi*m/(MI*dx);
k_Y_Rectangular=2*pi*n/(NI*dy);

Index = 1;
for f_Index = 100:1:110
    
    close all;
    f=freq(f_Index);

    for iy=1:1:N
       for ix=1:1:M 
           NF_X_Complex(ix,iy)=sdata.s21{iy,ix}(f_Index);
           NF_Y_Complex(ix,iy)=sdata2.s21{iy,ix}(f_Index);
       end
    end
    
    NF_X_Magnitude = 20*log10(abs(NF_X_Complex));
    NF_Y_Magnitude = 20*log10(abs(NF_Y_Complex));
    NF_Slots_X_Magnitude(Index,:) = interp2(x,y,NF_X_Magnitude',[0 0 0 0 0 0 0 0],0.0254*[-0.28 -0.20 -0.12 -0.04 0.04 0.12 0.20 0.28]);
    NF_Slots_Y_Magnitude(Index,:) = interp2(x,y,NF_Y_Magnitude',[0 0 0 0 0 0 0 0],0.0254*[-0.28 -0.20 -0.12 -0.04 0.04 0.12 0.20 0.28]);

    figure;
    subplot(2,1,1)
    surf(x*1000,y*1000,NF_X_Magnitude');
    title(sprintf('f = %f GHz',f/1000000000));
    xlabel('x (mm)');
    ylabel('y (mm)');
    zlabel('|E_{x}| (dB)');
    set(gca,'XLim',[min(x)*1000 max(x)*1000]);
    set(gca,'YLim',[min(y)*1000 max(y)*1000]);
    view(-37.5,30);
    shading flat;
    colorbar;
    subplot(2,1,2);
    surf(x*1000,y*1000,NF_Y_Magnitude');
    xlabel('x (mm)');
    ylabel('y (mm)');
    zlabel('|E_{y}| (dB)');
    set(gca,'XLim',[min(x)*1000 max(x)*1000]);
    set(gca,'YLim',[min(y)*1000 max(y)*1000]);
    view(-37.5,30);
    shading flat;
    colorbar;
    %print(gcf,'-dtiff',['NF_dB_' num2str(f) '_GHz']);
    clear NF_X_Magnitude NF_Y_Magnitude;

    [NF_X_Phase] = GoldsteinUnwrap2D(NF_X_Complex);
    [NF_Y_Phase] = GoldsteinUnwrap2D(NF_Y_Complex);
    NF_Slots_X_Phase(Index,:) = interp2(x,y,NF_X_Phase',0.0254*[-0.28 -0.20 -0.12 -0.04 0.04 0.12 0.20 0.28],[0 0 0 0 0 0 0 0]);
    NF_Slots_Y_Phase(Index,:) = interp2(x,y,NF_Y_Phase',0.0254*[-0.28 -0.20 -0.12 -0.04 0.04 0.12 0.20 0.28],[0 0 0 0 0 0 0 0]);
    
    figure;
    subplot(2,2,1)
    surf(x*1000,y*1000,NF_X_Phase');
    title(sprintf('f = %f GHz',f/1000000000));
    xlabel('x (mm)');
    ylabel('y (mm)');
    zlabel('\angle E_{x} (rad)');
    set(gca,'XLim',[min(x)*1000 max(x)*1000]);
    set(gca,'YLim',[min(y)*1000 max(y)*1000]);
    view(-37.5,30);
    shading flat;
    colorbar;
    subplot(2,2,2);
    imagesc(NF_X_Phase');
    colorbar;
    subplot(2,2,3);
    surf(x*1000,y*1000,NF_Y_Phase');
    xlabel('x (mm)');
    ylabel('y (mm)');
    zlabel('\angle E_{y} (rad)');
    set(gca,'XLim',[min(x)*1000 max(x)*1000]);
    set(gca,'YLim',[min(y)*1000 max(y)*1000]);
    view(-37.5,30);
    shading flat;
    colorbar;
    subplot(2,2,4);
    imagesc(NF_Y_Phase');
    colorbar;
    %print(gcf,'-dtiff',['NF_rad_' num2str(f) '_GHz']);
    clear NF_X_Phase NF_Y_Phase;

    % See equations (16-7a) and (16-7b) in Balanis
    f_X_Rectangular=ifftshift(ifft2(NF_X_Complex,MI,NI));
    f_Y_Rectangular=ifftshift(ifft2(NF_Y_Complex,MI,NI));
    
    f_X_Rectangular_Magnitude=20*log10(abs(f_X_Rectangular));
    f_Y_Rectangular_Magnitude=20*log10(abs(f_Y_Rectangular));

    figure;
    subplot(2,1,1);
    surf(k_X_Rectangular,k_Y_Rectangular,f_X_Rectangular_Magnitude');
    title(sprintf('f = %f GHz',f/1000000000));
    xlabel(sprintf('k_{x} (m^{-1})'));
    ylabel(sprintf('k_{y} (m^{-1})'));
    zlabel('|f_{x}| (dB)');
    set(gca,'XLim',[min(k_X_Rectangular) max(k_X_Rectangular)]);
    set(gca,'YLim',[min(k_Y_Rectangular) max(k_Y_Rectangular)]);
    view(-37.5,30);
    shading flat;
    colorbar;
    subplot(2,1,2);
    surf(k_X_Rectangular,k_Y_Rectangular,f_Y_Rectangular_Magnitude');
    xlabel(sprintf('k_{x} (m^{-1})'));
    ylabel(sprintf('k_{y} (m^{-1})'));
    zlabel('|f_{y}| (dB)');
    set(gca,'XLim',[min(k_X_Rectangular) max(k_X_Rectangular)]);
    set(gca,'YLim',[min(k_Y_Rectangular) max(k_Y_Rectangular)]);
    view(-37.5,30);
    shading flat;
    colorbar;

    lambda0=c/f;
    k0=2*pi/lambda0;
    f_X_Spherical=interp2(k_X_Rectangular,k_Y_Rectangular,abs(f_X_Rectangular'),k0*sin(theta).*cos(phi),k0*sin(theta).*sin(phi));
    f_Y_Spherical=interp2(k_X_Rectangular,k_Y_Rectangular,abs(f_Y_Rectangular'),k0*sin(theta).*cos(phi),k0*sin(theta).*sin(phi));

    r=10000;
    C=j*(k0*exp(-j*k0*r))/(2*pi*r);
    Etheta=C*(f_X_Spherical.*cos(phi)+f_Y_Spherical.*sin(phi));
    Ephi=C*cos(theta).*(-f_X_Spherical.*sin(phi)+f_Y_Spherical.*cos(phi));
    W=1/(2*120*pi).*(Etheta.*conj(Etheta)+Ephi.*conj(Ephi));

    figure;
    subplot(2,1,1);
    surf(theta,phi,20*log10(abs(f_X_Spherical)));
    title(sprintf('f = %f GHz',f/1000000000));
    xlabel('\theta (rad)');
    ylabel('\phi (rad)');
    zlabel('|f_{x}| (dB)');
    axis([-pi/2 pi/2 0 pi min(min(20*log10(abs(f_X_Spherical)))) max(max(20*log10(abs(f_X_Spherical))))]);
    view(-37.5,30);
    shading flat;
    colorbar;
    subplot(2,1,2);
    surf(theta,phi,20*log10(abs(f_Y_Spherical)));
    xlabel('\theta (rad)');
    ylabel('\phi (rad)');
    zlabel('|f_{y}| (dB)');
    axis([-pi/2 pi/2 0 pi min(min(20*log10(abs(f_Y_Spherical)))) max(max(20*log10(abs(f_Y_Spherical))))]);
    view(-37.5,30);
    shading flat;
    colorbar;

    figure;
    subplot(2,1,1);
    surf(theta,phi,20*log10(abs(Etheta)));
    title(sprintf('f = %f GHz',f/1000000000));
    xlabel('\theta (rad)');
    ylabel('\phi (rad)');
    zlabel('|E_{\theta}| (dBi)');
    axis([-pi/2 pi/2 0 pi min(min(20*log10(abs(Etheta)))) max(max(20*log10(abs(Etheta))))]);
    view(-37.5,30);
    shading flat;
    colorbar;
    subplot(2,1,2);
    surf(theta,phi,20*log10(abs(Ephi)));
    xlabel('\theta (rad)');
    ylabel('\phi (rad)');
    zlabel('|E_{\phi}| (dBi)');
    axis([-pi/2 pi/2 0 pi min(min(20*log10(abs(Ephi)))) max(max(20*log10(abs(Ephi))))]);
    view(-37.5,30);
    shading flat;
    colorbar;

    figure;
    subplot(1,2,1);
    sphere3d(20*log10(abs(Etheta))-max(max(20*log10(abs(Etheta')))),0,pi,-pi/2,pi/2,1,1,'surf','spline');
    title('|E_{\theta}| (dBi)');
    subplot(1,2,2);
    sphere3d(20*log10(abs(Ephi))-max(max(20*log10(abs(Ephi')))),0,pi,-pi/2,pi/2,1,1,'surf','spline');
    title('|E_{\phi}| (dBi)');

    figure;
    surf(theta,phi,10*log10(W));
    title(sprintf('f = %f GHz',f/1000000000));
    xlabel('\theta (rad)');
    ylabel('\phi (rad)');
    zlabel('|W| (dBi)');
    axis([-pi/2 pi/2 0 pi min(min(10*log10(W))) max(max(10*log10(W)))]);
    view(-37.5,30);
    shading flat;
    colorbar;

    W_Size=size(W);
    EPLANE(Index,:)=10*log10(W(floor(W_Size(2)/2),:))-max(10*log10(W(floor(W_Size(2)/2),:)));
    HPLANE(Index,:)=10*log10(W(1,:))-max(10*log10(W(1,:)));

    Index=Index+1;

end

figure;
subplot(1,2,1);
plot(180/pi*theta(1,:),EPLANE,[-25 -25],[-30 0],'k',[25 25],[-30 0],'k',[-90 90],[-13.6 -13.6],'k');
title('E-Plane');
xlabel('LOAD \leftarrow \theta (Deg) \rightarrow INPUT');
ylabel('Directivity (dBi)');
set(gca,'XLim',[-90 90]);
set(gca,'YLim',[-30 0]);
subplot(1,2,2);
plot(180/pi*theta(1,:),HPLANE);
title('H-Plane');
xlabel('\theta (Deg)');
ylabel('Directivity (dBi)');
set(gca,'XLim',[-90 90]);
set(gca,'YLim',[-30 0]);

save RESULTS EPLANE HPLANE NF_Slots_X_Magnitude NF_Slots_Y_Magnitude NF_Slots_X_Phase NF_Slots_Y_Phase freq