Commit 4e7f304d authored by Marco Berzborn's avatar Marco Berzborn

cleanup test files

parent cd96dc67
ccx
% <ITA-Toolbox>
% 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.
% </ITA-Toolbox>
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);
This diff is collapsed.
%% BOSCH Topic - Fourpole (two-port) in matrix notation
% Author: Pascal Dietrich - 2011 - pdi@akustik.rwth-aachen.d
% CONFIDENTIAL
%
% <ITA-Toolbox>
% 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.
% </ITA-Toolbox>
%% 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'})
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment