ccx % % This file is part of the application TPA-TPS for the ITA-Toolbox. All rights reserved. % You can find the license for this m-file in the application folder. % a = ita_generate('impulse',1,44100,16); b = a* 0; a_mat = [a, b; b*2, -2*a]; b_mat = [1*a, 1*a; 1*a, 1*a]; c = [1 1; 1 1]; %% Y matrix fft_degree = 16; folder = '/Users/pascaldietrich/MATLAB/__Bosch/BOSCH - auralizationBox otpa/vonMLI/Y_r'; for idx = 1:3 for jdx = 1:3 Yr(idx,jdx) = ita_extract_dat(ita_read([folder filesep 'y' num2str(idx) num2str(jdx) '.ita']),fft_degree,'symmetric'); end end folder = '/Users/pascaldietrich/MATLAB/__Bosch/BOSCH - auralizationBox otpa/vonMLI/Ys_model'; for idx = 1:3 for jdx = 1:3 Ys(idx,jdx) = ita_extract_dat(ita_mpb_filter( ita_read([folder filesep 'y' num2str(idx) num2str(jdx) '.ita']),[10 0], 'zerophase'),fft_degree,'symmetric'); end end %% fit tic for idx = 1:3 for jdx = 1:3 x = ita_frequency_dependent_time_window(Ys(idx,jdx),[0.6 0.7; 0.2 0.3; 0.05 0.06], [1000 3000], 'symmetric'); Ys_fit(idx,jdx) = ita_audio2zpk_rationalfit(x,'degree',80,'freqRange',[30 8000]); end end for idx = 1:3 for jdx = 1:3 x = ita_frequency_dependent_time_window(Yr(idx,jdx),[0.6 0.7; 0.2 0.3; 0.05 0.06], [1000 3000], 'symmetric'); Yr_fit(idx,jdx) = ita_audio2zpk_rationalfit(Yr(idx,jdx),'degree',80,'freqRange',[30 8000],'mode','log'); end end toc %% write folder = 'E:\pdi_daten\MATLAB\TPA-TPS'; ita_write(Ys_fit,[folder filesep 'Ys_fit.ita']) ita_write(Yr_fit,[folder filesep 'Yr_fit.ita']) %% load folder = 'E:\pdi_daten\MATLAB\TPA-TPS'; Ys_fit = load([folder filesep 'Ys_fit.ita'],'-mat'); Yr_fit = ita_read([folder filesep 'Yr_fit.ita']); %% for idx = 1:3 for jdx = 1:3 Ys_test(idx,jdx) = Ys_fit(idx,jdx)';% / Ys(idx,jdx); Yr_test(idx,jdx) = Yr_fit(idx,jdx)';% / Yr(idx,jdx); end end %% coupling tic K = Ys / (Yr + Ys); toc %% % K_fit = Ys_fit / (Yr_fit + Ys_fit); for idx = 1:3 for jdx = 1:3 % K_fit(idx,jdx).channelNames{1} = ['K' num2str(idx) num2str(jdx)]; K(idx,jdx).channelNames{1} = ['K' num2str(idx) num2str(jdx)]; end end %% K_test = Ys_test / (Yr_test + Ys_test); %% for idx = 1:3 for jdx = 1:3 K(idx,jdx).channelNames{1} = ['K' num2str(idx) num2str(jdx)]; end end ita_plot_spkphase(merge(K)) %% res1 = a_mat * c; ita_plot_spk(merge(res1),'nodb','ylim',[-5 5]) res2 = a_mat * b_mat; ita_plot_spk(merge(res2),'nodb','ylim',[-5 5]) %% coupling Ys = [1 0; 0, 1]; Yr = [0.9 0.1; 0.1 0.9]; Ys_mat = [1*a, 0*a; 0*a, 1*a]; Yr_mat = [.9*a, .1*a; .1*a, .9*a]; K = Ys / (Yr + Ys) for idx = 1:3 for jdx = 1:3 Yr(idx,jdx) = ita_extract_dat(ita_read([folder filesep 'y' num2str(idx) num2str(jdx) '.ita']),fft_degree,'symmetric'); end end K_mat = Ys_mat / (Yr_mat + Ys_mat) %% ccx %% mli raw TP with hammer clear TP for foot_idx = 1:3 for idx = [1] data = ita_read(['\\verdi\scratch\lievens\pdi\steel_cones\tps\Yc_foot' num2str(foot_idx) '_MDF_25e-3__pvc_run' num2str(idx) '.ita']); F = ita_time_window(data.ch(1),[0.02 0 0.5 0.7],'time'); p = data.ch(2:6); TP = p/F; end res(foot_idx) = ita_time_window(TP,[0.9 1.1],'time'); end %% mli data = ita_read('Z:\lievens\pdi\steel_cones\Fvp__4-20Hz_expsweep_fft20_8WashT.ita'); a = data.ch([2 3 1]); p = data.ch([7:11]); F = data.ch(4:6); % imp = F.ch(1)*0; % imp.timeData(:,1) = imp.nSamples; % F = merge(F,imp); %% OPA TP = ita_otpa(p,F,'blocksize',4096*8,'overlap',0.5,'tol',0.01,'window',@hann); TPopa = merge(TP.ch(1)); TPopa.plot_spk x = ita_time_window(TPopa,[0.002 0 double(TPopa.trackLength)*[0.9 0.99]],'time'); x.plot_spk %% pre-white spectrum TP = ita_otpa(p,F,'blocksize',4096*8,'overlap',0.5,'tol',0.1,'window',@hann,'prewhite'); for idx = 1:TP(1).nChannels TPopa(idx) = merge(TP.ch(idx)); end % TPopa(2).plot_spk % x = ita_time_window(TPopa,[0.002 0 double(TPopa.trackLength)*[0.5 0.99]],'time'); % x.plot_spk %% comparison close all idx = 1; TP(idx).plot_spk('ylim',[-40 0]); x = merge(res.ch(idx)); x.plot_spk('ylim',[-40 0]); %% synthesis for pidx = 1:numel(TP) p_test(pidx) = ita_sum(F * ita_extend_dat(TP(pidx),a.nSamples,'symmetric')); end p_test = merge(p_test); %% hammer messung data_hammer = ita_read('Z:\lievens\pdi\steel_cones\Yc_foot1_MDF_25e-3__pvc_run1.ita'); data_hammer = ita_time_shift(data_hammer,'30dB'); data_hammer = ita_time_shift(data_hammer,0.025,'time'); data_hammer = ita_time_window(data_hammer,[0.02 0],'time'); data_hammer = ita_frequency_dependent_time_window(data_hammer,[0.5 1; 0.05 0.1],500); res = ita_divide_spk( data_hammer.ch([1:3 5:12]),data_hammer.ch(4),'regularization',[20 4000]); TPorig = ita_divide_spk(res,res.ch(4),'regularization',[20 4000]); TPorigF = TPorig.ch(7:11); TPorigF.plot_spk %% synthesis with orig for pidx = 1:numel(TP) p_test_orig(pidx) = ita_sum(F * ita_extend_dat(TPorigF.ch(pidx),a.nSamples,'symmetric')); end p_test_orig = merge(p_test_orig); %% OPA with synthesis data TP = ita_otpa(p_test_orig,F,'blocksize',4096,'overlap',0.5,'tol',0.00000000001,'window',@hann); TPopa = merge(TP.ch(1)) TPopa.plot_spk %% OPA with randomized phase TPrand = ita_otpa(ita_randomize_phase(p),ita_randomize_phase(F),'blocksize',4096*8,'overlap',0.5,'tol',0.00000000001,'window',@hann); TPoparand = merge(TPrand.ch(1)) TPoparand.plot_spk %% ************************************************************************ %% Excitation Signal F0 = itaAudio(); F0.samplingRate = 100; F0.time = zeros(100,1); F0.time = 0.2*sin(2*F0.timeVector* pi); F0.time(1) = 1; F0.time(20) = -1; n = itaAudio; n.trackLength = 15; n.time = (n.timeVector / n.trackLength)*50 + 20; % n.time = 500 * n.time * 0; F1 = F0; F1.time = sin(1*F1.timeVector*2*pi) + 0.3 * sin(2*F1.timeVector*2*pi) + cos(4*F1.timeVector*2*pi); n1 = itaAudio; n1.trackLength = double(n.trackLength); n1.time = (n.timeVector / n1.trackLength)*1300 + 100+ 20*sin(40*n.timeVector/double(n1.trackLength) * 2 * pi); excitationsignal(1) = ita_normalize_dat( ita_iem_force_transform (F0, n,'oversample',100,'periodic',true)); excitationsignal(2) = 0.1 * ita_normalize_dat( ita_iem_force_transform (F1, n1)); %% amp = 1; sr = n.samplingRate; nSamples = n1.samplingRate * double(n1.trackLength); freq_vec = [20 10000]; stop_margin = 1; NTIcoeffs = [1 0.2]; % Generating a long MLS signal mls_raw = ita_generate('mls', amp, sr, 20); % Performing a hole in the MLS signal mls = mls_raw; mls.timeData(1:30000) = 0; % Sweep and compensation sweep = amp*ita_generate('linsweep',freq_vec,0.1,sr,nSamples); sweep = ita_extend_dat(sweep,mls.nSamples); % sweep = ita_time_shift(sweep,stop_margin/4); sweep.signalType = 'energy'; % because acts as a filter silencesweep = mls*sweep; silencesweep = ita_extend_dat(silencesweep,nSamples); % Silence Sweep excitationsignal(3) = ita_normalize_dat(silencesweep); %% expsweep %% playback test = merge(excitationsignal); %% ita_portaudio(test,'OutputChannels',[3 4 1]) %% ita_portaudio(test.ch(3),'OutputChannels',[1]) %% MS MS = ita_measurement %% MS messung ex = ita_amplify( MS.excitation ,'-15dB'); ex.trackLength = 10; ex = merge(ex,0.5*ita_time_shift(ex,3,'time'),ita_time_shift(ex,6,'time')); ex = ita_time_shift(ex,1,'time'); %% ita_portaudio(ex,'OutputChannels',[3 4 1]); %% ita_portaudio(ex.ch(1),'OutputChannels',[1]); %% test motor out_ch = [3 2]; motor_sweep = ita_generate('expsweep',[5000 17000],0.1,sr,nSamples); ita_portaudio(motor_sweep,'OutputChannels',out_ch) %% folder folder = '/Users/pascaldietrich/MATLAB/BOSCH - auralizationBox otpa'; cd (folder) %% auswertung F_ch = 1:3; p_ch = 13:15; TP = ita_otpa(data.ch(p_ch),data.ch(F_ch),'blocksize',4096,'overlap',0.5,'tol',0.00000000001,'window',@hann); %% deconv ir = ita_divide_spk(data,data.ch(14),'regularization',[20 10000]); ir_win = ita_time_window(ir,[0.5 1],'symmetric'); %% ls sweep hexaeder ls_sweep = ita_generate('expsweep',[20 17000],0.1,19); ita_portaudio( ita_amplify( ls_sweep ,'-20dB'), 'OutputChannels', 3);