% Octave用   2005/11/25
%
% CSVファイルのデータを読取り,周波数分析し伝達関数を求める．
%

clear
format long e

flag_sampling=1;    %サンプリング周期の自動指定. '1':自動．'0':手動

col_u=2;            %入力 (分母)列の指定
col_y=3;            %出力 (分子)列の指定　
msec2sec=1000;     %時間が[msec]なら'1000'を[sec]なら'1'を代入．

filename00='D:\\05kadai_03\\98Test\\05_11\\05_11_24NC_real_test3\\';
filename01='051124_01g_No2freq.CSV';
loadfile=strcat(filename00,filename01)

testdata = dlmread(loadfile,',');	%周波数分析するデータファイル

start_number = 1;
end_number = size(testdata,1);

time_testdata = testdata(1:end_number,1);

data1 = testdata(1:(end_number-start_number+1),col_u);
data2 = testdata(1:(end_number-start_number+1),col_y);

%サンプリング周期[sec]の設定
	if (flag_sampling==1)
    		line01=3;
    		line02=103;
    		ij=1;
    		for ii=line01:line02
        		err(ij,1)=time_testdata(ii+1,1)-time_testdata(ii,1);
        		ij=ij+1;
    		end
    		TT = sum(err)/length(time_testdata(line01:line02,1))/msec2sec;
	elseif(flag_sampling==0)
   		TT = 100e-6;           %サンプリング周期[sec]
	end
	F=1/TT;     		         %サンプリング周波数[Hz]

MAG_max= 20;    %最大値(出力用)
MAG_min=-40;        %最小値(出力用)  

frq_max=100;%F/2; %最大周波数(出力用)
frq_min=0;        %最小周波数(出力用)

y01=fft(data1);
y02=fft(data2);

[Row,Col]=size(data1);

MAG01=abs(y01)/(Row/2);
MAG02=abs(y02)/(Row/2);


PHASE01= angle(y01)*180/pi;
PHASE02= angle(y02)*180/pi;


 Gain_dB = [];
 Phase_deg = [];
for ii = 1:Row

 Gain_dB = [Gain_dB
            20*log10(MAG02(ii)/MAG01(ii))];

 Phase_deg_tmp = PHASE02(ii)-PHASE01(ii);

	if( abs(Phase_deg_tmp) < 180 )
	Phase_deg = [Phase_deg
        	     Phase_deg_tmp ];
	elseif( Phase_deg_tmp < 0 )
 	Phase_deg = [Phase_deg
        	      (180-mod(-Phase_deg_tmp,180)) ];
	elseif( Phase_deg_tmp > 0 )
 	Phase_deg = [Phase_deg
        	      (mod(Phase_deg_tmp,180)-180) ];
	end

end

freq = (0:length(y01)-1)'*F/length(y01);

subplot(2,1,1)
clearplot();
subplot(2,1,2)
clearplot();

gset('lmargin 13')

subplot(2,1,1)
%axis([frq_min frq_max MAG_min MAG_max]);
axis([1 100 -40 20]);
xlabel('Frequency[Hz]')
ylabel('Gain[dB]')
grid('on')

plot(freq,Gain_dB,'-;;')

subplot(2,1,2)
xlabel('freqency[Hz]')
ylabel('Phase [deg]')
axis([1 100 -180 180] );
grid('on')

plot(freq,Phase_deg,'-;;')

save_data = [freq Gain_dB Phase_deg];

save_filename='051124_01g_No2freq_reslt.CSV';
save_filename=strcat(filename00,save_filename)

dlmwrite(save_filename,save_data,',');