clear;
clc;

fprintf('Bode diagram of the Fourier-box\n\n');
fprintf('Figure 1: Bodediagram average measurement+SD\n\n');
fprintf('N.B. SD=0 is not displayed in the graph!\n\n');

fprintf('\nWARNING!This script looks for the file in the working directory.\n\n');

f_name=input('Filename: ','s');

[met]=dlmread(f_name,' ',0,0);

if mean(met(1,:))==0;
    nr_met=0;
end;

if mean(met(1,:))~=0;
    size_met=size(met);
    nr_met=(size_met(1)-1)/2;
    nr_punt=size_met(2);
    f_m=met(1,:);
    
    for i=1:nr_met;
        for j=1:nr_punt;
            amp_met(i,j)=met(i+1,j);
            fase_met(i,j)=met(i+1+nr_met,j);
        end;
    end;
    if nr_met~=1;
        for i=1:nr_punt;
            av_amp(i)=mean(amp_met(:,i));
            av_log_amp(i)=20*log10(av_amp(i));
            av_fase(i)=mean(fase_met(:,i));
        end;
        sd_amp=std(amp_met,0,1);
        for i=1:nr_punt;
            sd_log_amp(i)=2/log(av_amp(i))*sd_amp(i);
        end;
        sd_fase=std(fase_met,0,1);
    else 
        av_log_amp=20*log10(amp_met);
        av_fase=fase_met;
        sd_log_amp=zeros(1,nr_punt);
        sd_fase=zeros(1,nr_punt);
    end;
end;

min_amp=floor(min(av_log_amp-abs(sd_log_amp)));
max_amp=ceil(max(av_log_amp+abs(sd_log_amp)));
min_fase=floor(min(av_fase-abs(sd_fase)));
max_fase=ceil(max(av_fase+abs(sd_fase)));

if nr_met~=0;
    figure(1);
    set(clf,'PaperType','A4');
    set(clf,'Color',[1,1,1]);
    subplot(2,1,1);
    graf=semilogx(f_m,av_log_amp,'.-b','LineWidth',1);
%    legend([graf],'meting',3);
    hold on;
    width_sd=(log10(f_m(2))-log10(f_m(1)))/3;
    for i=1:nr_punt;
        if sd_log_amp(i)~=0;
            x=[f_m(i),f_m(i)];
            y=[av_log_amp(i)+sd_log_amp(i),av_log_amp(i)-sd_log_amp(i)];
            xl1=[10^(log(f_m(i))/log(10)-width_sd),10^(log(f_m(i))/log(10)+width_sd)];
            yl1=[av_log_amp(i)+sd_log_amp(i),av_log_amp(i)+sd_log_amp(i)];
            yl2=[av_log_amp(i)-sd_log_amp(i),av_log_amp(i)-sd_log_amp(i)];
            plot(x,y,'-b');
            plot(xl1,yl1,'-b');
            plot(xl1,yl2,'-b');        
        end;
    end;
    hold off;
    ylabel('Amplitude [dB]','Fontsize',10);

    % grid bepaling
    startfm=floor(log10(f_m(1)));
    sizefm=size(f_m);
    endfm=ceil(log10(f_m(sizefm(2))));
%    xfm(1)=10^startfm;
    xfm(1)=f_m(1);
    i=startfm+1;
    countm=2;
    flag=0;
    while flag~=1;
        if xfm(countm-1)<f_m(nr_punt);
            xfm(countm)=10^(i);
            countm=countm+1;
            i=i+1;
        else flag=1;
        end;
    end;
    xfm(countm-1)=f_m(nr_punt);
%    set(gca,'XTick',xfm,'Fontsize',9);
    
    min_amp_met=floor(min(av_log_amp-abs(sd_log_amp)))-2;
    max_amp_met=ceil(max(av_log_amp+abs(sd_log_amp)))+2;
    
    set(gca,'Fontsize',9);
    axis([xfm(1),xfm(countm-1),min_amp_met,max_amp_met]);
    grid on;
    title('Bodediagram Fourierbox (Measurement)','Fontsize',10,'Fontweight','bold');

    subplot(2,1,2);
    
    graf=semilogx(f_m,av_fase,'.-b','LineWidth',1);
%    legend([graf],'meting',3);
    hold on;
    width_sd=(log10(f_m(2))-log10(f_m(1)))/3;
    for i=1:nr_punt;
        if sd_fase(i)~=0;
            x=[f_m(i),f_m(i)];
            y=[av_fase(i)+sd_fase(i),av_fase(i)-sd_fase(i)];
            xl1=[10^(log10(f_m(i))-width_sd),10^(log10(f_m(i))+width_sd)];
            yl1=[av_fase(i)+sd_fase(i),av_fase(i)+sd_fase(i)];
            yl2=[av_fase(i)-sd_fase(i),av_fase(i)-sd_fase(i)];
            plot(x,y,'-b');
            plot(xl1,yl1,'-b');
            plot(xl1,yl2,'-b');        
        end;
    end;
    hold off;
    
    xlabel('Frequency [Hz]','Fontsize',10);
    ylabel('Fase [rad]','Fontsize',10);
%    set(gca, 'XTick',xfm,'Fontsize',9);
 
    min_fase_met=floor(min(av_fase-abs(sd_fase)));
    max_fase_met=ceil(max(av_fase+abs(sd_fase)))+2;
 
    set(gca,'Fontsize',9);
    axis([f_m(1),f_m(nr_punt),min_fase,max_fase]);
    grid on;
end;