From 4e7f304d9e532e1c435f6927ca75fb88728a4887 Mon Sep 17 00:00:00 2001 From: Marco Berzborn Date: Thu, 17 Nov 2016 11:15:09 +0100 Subject: [PATCH] cleanup test files --- applications/TPA-TPS/test_pdi_otpa_tps.m | 332 ------------ .../TPA-TPS/test_pdi_plot_BOSCH_report.m | 480 ------------------ .../test_pdi_tpa_fourpole_matrix_notation.m | 169 ------ 3 files changed, 981 deletions(-) delete mode 100644 applications/TPA-TPS/test_pdi_otpa_tps.m delete mode 100644 applications/TPA-TPS/test_pdi_plot_BOSCH_report.m delete mode 100644 applications/TPA-TPS/test_pdi_tpa_fourpole_matrix_notation.m diff --git a/applications/TPA-TPS/test_pdi_otpa_tps.m b/applications/TPA-TPS/test_pdi_otpa_tps.m deleted file mode 100644 index 064c0314..00000000 --- a/applications/TPA-TPS/test_pdi_otpa_tps.m +++ /dev/null @@ -1,332 +0,0 @@ -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); - - - diff --git a/applications/TPA-TPS/test_pdi_plot_BOSCH_report.m b/applications/TPA-TPS/test_pdi_plot_BOSCH_report.m deleted file mode 100644 index 21aa42d4..00000000 --- a/applications/TPA-TPS/test_pdi_plot_BOSCH_report.m +++ /dev/null @@ -1,480 +0,0 @@ -%% *********************************************************************** -% ************************************************************************ -% ************************* plots , additional script ******************** -% ************************************************************************ -% ************************************************************************ -% -% Author Pascal Dietrich 2011 - pdi@akustik.rwth-aachen.de -% For BOSCH Report -% --- CONFIDENTIAL --- -% -% ************************************************************************ -% ************************************************************************ -% ************************************************************************ -% ************************************************************************ - -% -% 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. -% - - -%% Yc test -Yr_in_raw = ita_AudioAnalyticRationalMatrix2itaAudioMatrix(Yr_backup, 'fftDegree',fftDegree,'samplingRate',sr); -Ys_in_raw = ita_AudioAnalyticRationalMatrix2itaAudioMatrix(Ys_backup, 'fftDegree',fftDegree,'samplingRate',sr); -Yr_in_raw(1,1).comment = 'Yr'; -Ys_in_raw(1,1).comment = 'Ys'; -ita_disp('conversion done.') -ext = ''; -ylimits = [-60 -30]; - -%% all DOF -close all -comment = 'All DOFs'; -filename = comment(isstrprop(comment,'alpha')); -Yr2 = Yr_in_raw; -Ys2 = Ys_in_raw; -Yc = Yr2*pinv(Yr2 + Ys2)*Ys2; -Yc(1,1).comment = 'Yc'; -ita_tpa_plot_matrix({Ys2, Yr2,Yc},'freqVector',200,'filename',filename,'ticksF',{'FZ','MX','MY'} ,'ticksV',{'UZ','RX','RY'} ) -count = 0; -clear aux -for idx = 1:3:size(Yc,1) - count = count + 1; - aux(count) = Yc(idx,idx); -end - -Yc_noFoot_selected = aux.merge; -Yc_noFoot_selected.plot_spkphase('ylim',ylimits); -title(comment) -ita_savethisplot_gle('filename',[filename ext]); - -%cond number -cYr2 = ita_tpa_plot_matrix_condition_number(Yr2); -cYr2.channelNames = {'Yr'}; - -cYs2 = ita_tpa_plot_matrix_condition_number(Ys2); -cYs2.channelNames = {'Ys'}; - -cYc2 = ita_tpa_plot_matrix_condition_number(Yc); -cYc2.channelNames = {'Yc'}; - -ita_plot_spk(merge(cYs2, cYr2, cYc2),'nodb') - -%% all DOF - no foot coupling Yr Ys -close all -comment = 'All DOFs -- No Coupling Between Feet (Yr, Ys)'; -filename = comment(isstrprop(comment,'alpha')); -Yr2 = diag(Yr_in_raw,[],3); -Ys2 = diag(Ys_in_raw,[],3); -Yc = Yr2*pinv(Yr2 + Ys2)*Ys2; -Yc(1,1).comment = 'Yc'; -ita_tpa_plot_matrix({Ys2, Yr2,Yc},'freqVector',200,'filename',filename,'ticksF',{'FZ','MX','MY'} ,'ticksV',{'UZ','RX','RY'} ) -count = 0; -clear aux -for idx = 1:3:size(Yc,1) - count = count + 1; - aux(count) = Yc(idx,idx); -end - -Yc_noFoot_selected = aux.merge; -Yc_noFoot_selected.plot_spkphase('ylim',ylimits); -title(comment) -ita_savethisplot_gle('filename',[filename ext]); - -% all DOF - no foot coupling Yr -close all -comment = 'All DOFs -- No Coupling Between Feet (Yr)'; -filename = comment(isstrprop(comment,'alpha')); -Yr2 = diag(Yr_in_raw,[],3); -Ys2 = Ys_in_raw; -Yc = Yr2*pinv(Yr2 + Ys2)*Ys2; -Yc(1,1).comment = 'Yc'; -ita_tpa_plot_matrix({Ys2, Yr2,Yc},'freqVector',200,'filename',filename,'ticksF',{'FZ','MX','MY'} ,'ticksV',{'UZ','RX','RY'} ) -count = 0; -clear aux -for idx = 1:3:size(Yc,1) - count = count + 1; - aux(count) = Yc(idx,idx); -end - -Yc_noFoot_selected = aux.merge; -Yc_noFoot_selected.plot_spkphase('ylim',ylimits); -title(comment) -ita_savethisplot_gle('filename',[filename ext]); - -% all DOF - no foot coupling Ys -close all -comment = 'All DOFs -- No Coupling Between Feet (Ys)'; -filename = comment(isstrprop(comment,'alpha')); -Yr2 = Yr_in_raw; -Ys2 = diag(Ys_in_raw,[],3); -Yc = Yr2*pinv(Yr2 + Ys2)*Ys2; -Yc(1,1).comment = 'Yc'; -ita_tpa_plot_matrix({Ys2, Yr2,Yc},'freqVector',200,'filename',filename,'ticksF',{'FZ','MX','MY'} ,'ticksV',{'UZ','RX','RY'} ) -count = 0; -clear aux -for idx = 1:3:size(Yc,1) - count = count + 1; - aux(count) = Yc(idx,idx); -end - -Yc_noFoot_selected = aux.merge; -Yc_noFoot_selected.plot_spkphase('ylim',ylimits); -title(comment) -ita_savethisplot_gle('filename',[filename ext]); - -% all DOF - no foot coupling - no DOF coupling -close all -comment = 'All DOFs -- No Coupling Between DOF'; -filename = comment(isstrprop(comment,'alpha')); -Yr2 = diag(Yr_in_raw,3); -Ys2 = diag(Ys_in_raw,3); -Yc = Yr2*pinv(Yr2 + Ys2)*Ys2; -Yc(1,1).comment = 'Yc'; -ita_tpa_plot_matrix({Ys2, Yr2,Yc},'freqVector',200,'filename',filename,'ticksF',{'FZ','MX','MY'} ,'ticksV',{'UZ','RX','RY'} ) -count = 0; -clear aux -for idx = 1:3:size(Yc,1) - count = count + 1; - aux(count) = Yc(idx,idx); -end -Yc_noFoot_selected = aux.merge; -Yc_noFoot_selected.plot_spkphase('ylim',ylimits); -title(comment) -ita_savethisplot_gle('filename',[filename ext]); - -% all DOF - no foot coupling - no DOF coupling -close all -comment = 'All DOFs -- No Coupling Between Feet (Yr, Ys) and Between DOF'; -filename = comment(isstrprop(comment,'alpha')); -Yr2 = diag(Yr_in_raw); -Ys2 = diag(Ys_in_raw); -Yc = Yr2*pinv(Yr2 + Ys2)*Ys2; -Yc(1,1).comment = 'Yc'; -ita_tpa_plot_matrix({Ys2, Yr2,Yc},'freqVector',200,'filename',filename,'ticksF',{'FZ','MX','MY'} ,'ticksV',{'UZ','RX','RY'} ) -count = 0; -clear aux -for idx = 1:3:size(Yc,1) - count = count + 1; - aux(count) = Yc(idx,idx); -end -Yc_noFoot_selected = aux.merge; -Yc_noFoot_selected.plot_spkphase('ylim',ylimits); -title(comment) -ita_savethisplot_gle('filename',[filename ext]); - -%% normal direction only - Should be the same as "all DOF - no foot coupling - no DOF coupling" -close all -comment = 'Normal Direction Only'; -filename = comment(isstrprop(comment,'alpha')); -Yr2 = diag(Yr_in_raw(1:3:end,1:3:end)); -Ys2 = diag(Ys_in_raw(1:3:end,1:3:end)); -Yc = Yr2*pinv(Yr2 + Ys2)*Ys2; -Yc(1,1).comment = 'Yc'; -ita_tpa_plot_matrix({Ys2, Yr2,Yc},'freqVector',200,'filename',filename,'ticksF',{'FZ'} ,'ticksV',{'UZ'} ) -count = 0; -clear aux -for idx = 1:1:size(Yc,1) - count = count + 1; - aux(count) = Yc(idx,idx); -end -Yc_noFoot_selected = aux.merge; -Yc_noFoot_selected.plot_spkphase('ylim',ylimits); -title(comment) -ita_savethisplot_gle('filename',[filename ext]); - -%% Ys invert due to ts -fmax = 2000; -sr = 5000; -close all -clc -ts = 5e-3; -Geoms = itaCoordinates([lxs lys ts]); %create itaCoordinates / dimensions -Ys = repmat(itaAudioAnalyticRational,coordFeet.nPoints*DOF,coordFeet.nPoints*DOF); -for idx = 1:coordFeet.nPoints - for jdx = 1:coordFeet.nPoints - Ys((idx-1)*DOF+(1:DOF),(jdx-1)*DOF+(1:DOF)) = ita_tps_mobility_plate(Geoms,coordFeet.n(jdx),coordFeet.n(idx),EX,PRXY,RHO,DMPR,'f_max',fmax); - end -end -nModes = numel(Ys(1,1).analyticData.f); -Ys(1,1) = 'Ys'; -Ys = Ys(idxF,idxF); -Ys = ita_AudioAnalyticRationalMatrix2itaAudioMatrix(Ys, 'fftDegree',fftDegree,'samplingRate',sr); -Ys.merge.plot_spkphase('ylim',[-40 70]) - - -fileStr = num2str(ts); -fileStr = fileStr(isstrprop(fileStr,'alphanum')); - -titleStr = ['ts: ' num2str(ts) ' - ' num2str(nModes) ' Modes']; -title(titleStr) -ita_savethisplot_gle('filename',['YSts' fileStr ]) - -ita_tpa_plot_matrix({Ys},'freqVector',200,'ticksF',{'FZ','MX','MY'} ,'ticksV',{'UZ','RX','RY'} ) -ita_plottools_aspectratio(1) -ita_savethisplot('matrix_0005.png') - - -ita_tpa_plot_matrix_condition_number(Ys); -title(titleStr) -ita_savethisplot_gle('filename',['condNumberts' fileStr ]) - -% % ita_tpa_plot_matrix({Ys},0,'freqVector',200,'ticksF',{'FZ','MX','MY'} ,'ticksV',{'UZ','RX','RY'} ) -% % ita_plottools_aspectratio(1) -% % print('-dpng','-r300',['matrix_' num2str(ts) '.png']) - - - -close all -clc -ts = 1e-3; -Geoms = itaCoordinates([lxs lys ts]); %create itaCoordinates / dimensions -Ys = repmat(itaAudioAnalyticRational,coordFeet.nPoints*DOF,coordFeet.nPoints*DOF); -for idx = 1:coordFeet.nPoints - for jdx = 1:coordFeet.nPoints - Ys((idx-1)*DOF+(1:DOF),(jdx-1)*DOF+(1:DOF)) = ita_tps_mobility_plate(Geoms,coordFeet.n(jdx),coordFeet.n(idx),EX,PRXY,RHO,DMPR,'f_max',fmax); - end -end -nModes = numel(Ys(1,1).analyticData.f); -Ys(1,1) = 'Ys'; -Ys = Ys(idxF,idxF); -Ys = ita_AudioAnalyticRationalMatrix2itaAudioMatrix(Ys, 'fftDegree',fftDegree,'samplingRate',sr); -Ys.merge.plot_spkphase('ylim',[-40 70]) - - -fileStr = num2str(ts); -fileStr = fileStr(isstrprop(fileStr,'alphanum')); - -titleStr = ['ts: ' num2str(ts) ' - ' num2str(nModes) ' Modes']; -title(titleStr) -ita_savethisplot_gle('filename',['YSts' fileStr ]) - -ita_tpa_plot_matrix_condition_number(Ys); -title(titleStr) -ita_savethisplot_gle('filename',['condNumberts' fileStr ]) - - - -close all -clc -ts = 0.1e-3; -Geoms = itaCoordinates([lxs lys ts]); %create itaCoordinates / dimensions -Ys = repmat(itaAudioAnalyticRational,coordFeet.nPoints*DOF,coordFeet.nPoints*DOF); -for idx = 1:coordFeet.nPoints - for jdx = 1:coordFeet.nPoints - Ys((idx-1)*DOF+(1:DOF),(jdx-1)*DOF+(1:DOF)) = ita_tps_mobility_plate(Geoms,coordFeet.n(jdx),coordFeet.n(idx),EX,PRXY,RHO,DMPR,'f_max',fmax); - end -end -nModes = numel(Ys(1,1).analyticData.f); -Ys(1,1) = 'Ys'; -Ys = Ys(idxF,idxF); -Ys = ita_AudioAnalyticRationalMatrix2itaAudioMatrix(Ys, 'fftDegree',fftDegree,'samplingRate',sr); -Ys.merge.plot_spkphase('ylim',[-40 70]) - - -fileStr = num2str(ts); -fileStr = fileStr(isstrprop(fileStr,'alphanum')); - -titleStr = ['ts: ' num2str(ts) ' - ' num2str(nModes) ' Modes']; -title(titleStr) -ita_savethisplot_gle('filename',['YSts' fileStr ]) - -ita_tpa_plot_matrix_condition_number(Ys); -title(titleStr) -ita_savethisplot_gle('filename',['condNumberts' fileStr ]) - - -%% condition number and damping -% close all -clc -sr = 5000; -fmax = 2000; -ts = 0.1e-3; -DMPR = 10; -% DMPR = 1e2; %damping - -Geoms = itaCoordinates([lxs lys ts]); %create itaCoordinates / dimensions -Ys = repmat(itaAudioAnalyticRational,coordFeet.nPoints*DOF,coordFeet.nPoints*DOF); -for idx = 1:coordFeet.nPoints - for jdx = 1:coordFeet.nPoints - Ys((idx-1)*DOF+(1:DOF),(jdx-1)*DOF+(1:DOF)) = ita_tps_mobility_plate(Geoms,coordFeet.n(jdx),coordFeet.n(idx),EX,PRXY,RHO,DMPR,'f_max',fmax); - end -end - -eigenFreq = Ys(1,1).analyticData.f; - -nModes = numel(Ys(1,1).analyticData.f); -Ys(1,1) = 'Ys'; -Ys = Ys(idxF,idxF); -Ys = ita_AudioAnalyticRationalMatrix2itaAudioMatrix(Ys, 'fftDegree',fftDegree,'samplingRate',sr); -% Ys.merge.plot_spkphase('ylim',[-40 70]) - -titleStr = ['ts: ' num2str(ts) ' - damping: ' num2str(DMPR) ' - ' num2str(nModes) ' Modes']; - -fileStr = [num2str(ts) '_' num2str(DMPR) ]; -fileStr = fileStr(isstrprop(fileStr,'alphanum')); - -ita_tpa_plot_matrix_condition_number(Ys,eigenFreq); -title(titleStr) -legend off -ylim([0 5e5]) - -ita_savethisplot_gle('filename',['condNumbertDamping' fileStr ]) - - -%% *********************** -%% *************** PLATE-BEAM ****************************************** -% % ccx -% This link has to be adapted! -% % load('F:\pdi_daten\MATLAB\BoschReport\ANSYS_plate_beam\workspace_data.mat') - -% unitsit -mulfac = 3; -v_units = [repmat(itaValue(1,'m/s'),1,mulfac) repmat(itaValue(1,'1/s'),1,mulfac)]; -F_units = [repmat(itaValue(1,'N'),1,mulfac) repmat(itaValue(1,'N m'),1,mulfac)]; -units = repmat(itaValue,2*mulfac,2*mulfac); % just init - -for idx = 1:numel(v_units) - for jdx = 1:numel(F_units) - aux = itaValue(v_units(idx)) / itaValue(F_units(jdx)); - units(jdx,idx) = aux; - end -end -units = repmat(units,3,3) - -%% get velocity -DOF = 6; -for idx = 1:DOF*3 - for jdx = 1:DOF*3 - Placa (idx,jdx) = ita_differentiate( Placa (idx,jdx)); - Placa (idx,jdx).channelUnits = units(idx,jdx); - PlacaViga (idx,jdx) = ita_differentiate( PlacaViga (idx,jdx)); - PlacaViga (idx,jdx).channelUnits = units(idx,jdx); - VigaMobility (idx,jdx) = ita_differentiate( VigaMobility (idx,jdx)); - VigaMobility (idx,jdx).channelUnits = units(idx,jdx); - end -end - -%% select DOF -selectDOF = 1:6; -% selectDOF = [3 4 5]; -idxF = ita_tpa_DOF2index(selectDOF,3,DOF); -idxR = ita_tpa_DOF2index(selectDOF,3,DOF); -Ys = VigaMobility(idxF,idxR); -Yr = Placa(idxF,idxR); -Yc = PlacaViga(idxF,idxR); -Yc_sim = Yr*pinv(Yr + Ys,1e-5)*Ys; -% Yc2_sim = Yr2*pinv(Yr2 + Ys2)*Ys2; - -%% naming -% Yc = Yc(idxF, idxR); -Yc(1,1).comment = 'Yc ANSYS'; -% Yc_sim = Yc_sim(idxF, idxR); -Yc_sim(1,1).comment = 'Yc calculated'; - -%% save -save(['admittance_data_plate_beam_' num2str(numel(selectDOF)) 'DOF.mat'],'Yr','Ys','Yc','Yc_sim') - -%% -if numel(selectDOF,3) - ita_tpa_plot_matrix({Ys, Yr, Yc, Yc_sim},'filename',['ansys_numerik_' num2str(numel(selectDOF)) 'DOF'],'ticksF',{'FZ','MX','MY'} ,'ticksV',{'UZ','RX','RY'},'freqVector',350) -else - ita_tpa_plot_matrix({Ys, Yr, Yc, Yc_sim},'filename',['ansys_numerik_' num2str(numel(selectDOF)) 'DOF'],'freqVector',350) -end - -%% windowing -Ys_win = ita_matrixfun(@ita_time_window,Ys,[0.4 0.6],'time','symmetric'); -Yr_win = ita_matrixfun(@ita_time_window,Yr,[0.4 0.6],'time','symmetric'); - -Yc_sim_win = Yr_win*pinv(Yr_win + Ys_win,1e-5)*Ys_win; - -%% -ita_tpa_plot_matrix_condition_number(Yc) -legend off -ylim([0 4e8]) -ita_savethisplot_gle('ansys_condnumber_yc') - -%% -ita_tpa_plot_matrix_condition_number(Ys) -legend off -ylimits = ylim; -ylim([0 1e11]) -ita_savethisplot_gle('ansys_condnumber_ys') - -%% -ita_tpa_plot_matrix_condition_number(Yr) -legend off -ylimits = ylim; -ylim([0 4e10]) -ita_savethisplot_gle('ansys_condnumber_yr') - -%% -count = 0 -clear aux -for idx = 1:3:size(Ys,1) - count = count + 1; - aux(count) = Yc_sim(idx,idx); -end -aux.merge.plot_spkphase(); - -%% Ys fit -for idx = 1:size(Ys,1) - for jdx = 1:size(Ys,2) - Ys_fit(idx,jdx) = ita_audio2zpk_rationalfit((Ys(idx,jdx)),'degree',15,'freqRange',[100 2000],'tendstozero',false); - end -end -Ys2 = ita_AudioAnalyticRationalMatrix2itaAudioMatrix(Ys_fit, 'fftDegree',Ys(1,1).nSamples,'samplingRate',4000); - -%% -for idx = 1:size(Ys,1) - for jdx = 1:size(Ys,2) - Ys(idx,jdx).signalType = 'energy'; - Ys_fit(idx,jdx).fftDegree = 16; - Yr_fit(idx,jdx).fftDegree = 16; - - Ys_fit(idx,jdx).samplingRate = Ys(1,1).samplingRate; - Yr_fit(idx,jdx).samplingRate = Ys(1,1).samplingRate; - - end -end - -%% -for idx = 1:size(Ys,1) - disp([num2str(idx) ' of ' num2str(size(Ys,1))]) - for jdx = 1:size(Ys,2) - Yr_fit(idx,jdx) = ita_audio2zpk_rationalfit((Yr(idx,jdx)),'degree',[200],'freqRange',[20 2000],'tendstozero',false); - end -end -Yr2 = ita_AudioAnalyticRationalMatrix2itaAudioMatrix(Yr_fit, 'fftDegree',Ys(1,1).nSamples,'samplingRate',4000); - -%% transpose? -Yr_trans = ita_tpa_matrix_transpose(Yr); -test = Yr.merge/Yr_trans.merge - -%% -data = Yr -ita_plot_spkphase(merge(data(1,1),data(4,4), data(7,7))) -legend off -ylim([-80 20]) -ita_savethisplot_gle('ansys_spkphase_Yr') - -%% -data = Yr_fit; -ita_plot_spkphase(merge(data(1,1)',data(4,4)', data(7,7)')) -legend off -ylim([-80 20]) -ita_savethisplot_gle('ansys_spkphase_Yr_fit') - -%% -Yr_fit2 = Yr_fit'; -Ys_fit2 = Ys_fit'; -Yc_fit = Yr_fit2 * pinv( Yr_fit2 + Ys_fit2) * Ys_fit2; - - - - - - - - diff --git a/applications/TPA-TPS/test_pdi_tpa_fourpole_matrix_notation.m b/applications/TPA-TPS/test_pdi_tpa_fourpole_matrix_notation.m deleted file mode 100644 index 9ad81fd3..00000000 --- a/applications/TPA-TPS/test_pdi_tpa_fourpole_matrix_notation.m +++ /dev/null @@ -1,169 +0,0 @@ -%% BOSCH Topic - Fourpole (two-port) in matrix notation -% Author: Pascal Dietrich - 2011 - pdi@akustik.rwth-aachen.d -% CONFIDENTIAL -% - -% -% 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. -% - - -%% Init -ccx -sr = 44100; -fftDegree = 14; - -%% Generate Scenario - Source: mass - Fourpole: additional mass - Receiver: spring -s = itaAudio(); -s.fftDegree = fftDegree; -s.samplingRate = sr; -s.freq = 1i * 2 * pi * s.freqVector; -s = s * itaValue('Hz'); - -%% Inits -Null = repmat(0*s/itaValue('Hz'),2,2); -Ys = Null; -Yr = Null; - -%% Source -m1 = itaValue('0.1 kg'); -m2 = itaValue('0.3 kg'); -Ys(1,1) = 1 / (s * m1); -Ys(2,2) = 1 / (s * m2); - -%% Receiver -n1 = itaValue(1e-7, 'm/N'); -n2 = itaValue(1e-8, 'm/N'); -Yr(1,1) = s * n1; -Yr(2,2) = s * n2; - -%% resonances without two-port -ita_disp('Resonance frequencies without two port') -f1 = 1/(2 * pi * sqrt( (m1) * n1)) -f2 = 1/(2 * pi * sqrt( (m2) * n2)) - -%% Yc without two-port -Yc = Yr * pinv(Yr + Ys) * Ys; -Yc(1,1).comment = 'Yc - no two-port'; -Yc.merge.plot_spkphase - -%% clean up -close all -clc - -%% Two-port with mass -delta_m1 = itaValue(0.01,'kg'); -delta_m2 = itaValue(0.01,'kg'); -A = 0*Null+eye(2); -B = Null * itaValue('kg/s'); -B(1,1) = - s * delta_m1; -B(2,2) = - s * delta_m2; -C = Null * itaValue('s/kg'); -D = 0*Null+eye(2); -T1 = [A B; C D]; - -%% Two Port with spring -delta_n1 = itaValue(1e-5, 'm/N'); -delta_n2 = itaValue(1e-5, 'm/N'); -A = 0*Null+eye(2); -B = Null * itaValue('kg/s'); -% B(1,1) = - 1/(s * delta_n1); -% B(2,2) = - 2/(s * delta_n2); -C = Null * itaValue('s/kg'); -C(1,1) = -(s * delta_n1); -C(2,2) = -(s * delta_n2); -D = 0*Null+eye(2); -T2 = [A B; C D]; - -%% Multiplication of single two-ports -T3 = T1 * T2 * T1; -A = T3(1:2,1:2); -B = T3(1:2,3:4); -C = T3(3:4,1:2); -D = T3(3:4,3:4); - -%% resonances including two-port -ita_disp('Resonance frequencies including two port') -f1 = 1/(2 * pi * sqrt( (m1 + delta_m1 ) * (n1 + delta_n1))) -f2 = 1/(2 * pi * sqrt( (m2 + delta_m2 ) * (n2 + delta_n2))) - -%% two-port transformation equation - now isolator is included into the source -Ys_incl = ita_tps_two_port_transformation(Ys, A, B, C, D); -Ys_incl.merge.plot_spkphase - -%% Ys_incl connected to structure Yr -Yc_tp = Yr * pinv(Yr + Ys_incl) * Ys_incl; -Yc_tp(1,1).comment = 'Yc - including two-port'; -Yc_tp.merge.plot_spkphase - -%% ********************** Isolators of AuralizationBox ******************** -%% Inits -s = itaAudio(); -s.fftDegree = fftDegree; -s.samplingRate = sr; -s.freq = 1i * 2 * pi * s.freqVector; -s = s * itaValue('Hz'); -Null = repmat(0*s/itaValue('Hz'),3,3); - -%% units INIT -F_units = {'N','N m','N m'}; -v_units = {'m/s','1/s','1/s'}; -for idx = 1:numel(F_units) - for jdx = 1:numel(v_units) - Aunit(idx,jdx) = itaValue(0,ita_deal_units(F_units{idx}, F_units{jdx},'/')); - Bunit(idx,jdx) = itaValue(0,ita_deal_units(F_units{idx}, v_units{jdx},'/')); - Cunit(idx,jdx) = itaValue(0,ita_deal_units(v_units{idx}, F_units{jdx},'/')); - Dunit(idx,jdx) = itaValue(0,ita_deal_units(v_units{idx}, v_units{jdx},'/')); - end -end -Ainit = Null * Aunit + eye(3); -Binit = Null * Bunit; -Cinit = Null * Cunit; -Dinit = Null * Dunit + eye(3); - -%% Two-port with mass / moment of inertia -% values -m = itaValue(0.01,'kg'); -r = itaValue(0.005,'m'); -l = itaValue(0.005,'m'); -J = m *(r^2/4 + l^2/12); - -A = Ainit;B = Binit; C=Cinit;D = Dinit; - -B(1,1) = - s * m; -B(2,2) = - s * J; -B(3,3) = - s * J; - -T_mass = [A B; C D]; - - -%% Two Port with spring -%values -n = itaValue(1e-5, 'm/N'); -n_rot = itaValue(1e-5, '1/N m'); - -A = Ainit;B = Binit; C=Cinit;D = Dinit; - -C(1,1) = -(s * n); -C(2,2) = -(s * n_rot); -C(3,3) = -(s * n_rot); - -T_spring = [A B; C D]; - -%% Isolator -T_iso = T_mass * T_spring * T_mass; -A_iso = T3(1:2,1:2); -B_iso = T3(1:2,3:4); -C_iso = T3(3:4,1:2); -D_iso = T3(3:4,3:4); - -m = itaValue(0.01,'kg'); -n = itaValue(1e-5, 'm/N'); -[A B C D] = ita_tpa_isolator(m,n,fftDegree,sr,{'N','N m','N m'},{'m/s','1/s','1/s'}); - -%% -m = itaValue(0,'kg'); -n = itaValue(0, 'm/N'); -[A B C D] = ita_tpa_isolator(m,n,fftDegree,sr,{'N','N m','N m'},{'m/s','1/s','1/s'}) - -- GitLab