Function inside function using lsqcurvefit

2 views (last 30 days)
Hello,
I want to fit the curve of the TL function to the xn_fft function and find the design variables x []. Notice that the TL function contains several functions. The code gives:
{Index exceeds matrix dimensions. Error in optimization (line 34). Ln1=x(1);Ln2=x(2); % necks length }
Any help is welcome and appreciated!
clear all;clc;
%% Harmonic sine waves
fs = 400; % Sampling rate [Hz]
Ts = 1/fs; % Sampling period [s]
fNy = fs / 2; % Nyquist frequency [Hz]
duration = 0.5; % Duration [s]
t = 0 : Ts : duration-Ts; % Time vector
noSamples = length(t); % Number of samples
f = 0 : fs/noSamples : fs - fs/noSamples; % Frequency vector
fundamental_signal = 1.5.*sin(2 .* pi .* 50 .* t); %Fundamental signal
First_Harmonic = 1.5.*sin(2 .* pi .* 100 .* t); %FH1=100 Hz %First Harmonic
Total_signal = fundamental_signal + First_Harmonic; % Contaminated signal
xn_fft = abs(fft(Total_signal)); % FFT
%% TL for 2-DOF HR
c=340.3; % Speed of sound in the air
w=2*pi()*f;
density=1.225; % Air density
Sd=4.3*4.3e-4; % duct area
Lv1=2.5e-3;Lv2=2.5e-3; %length of the cicular cavities
%Variables
x=[ ];
Ln1=x(1);Ln2=x(2); % necks length
S1=x(3); S2=x(4); % necks area
V1=x(5); V2=x(6); % Cavities volume
Rn1= (S1/(pi()))^0.5; Rn2= (S2/(pi()))^0.5; % necks area
Rv1=((V1/(pi()))^0.5)/0.05; Rv2=((V2/(pi()))^0.5)/0.05;
% Calculations
% Effective neck length
correc_n1d=0.82*(Rn1/2);
correc_n1v1=0.85*Rn1*(1-(1.25*(Rn1/Rv1)));
correc_n2v1=0.85*Rn2*(1-(1.25*(Rn2/Rv1)));
correc_n2v2=0.85*Rn2*(1-(1.25*(Rn2/Rv2)));
Leff1=Ln1+correc_n1d+correc_n1v1;
Leff2=Ln2+correc_n2v1+correc_n2v2;
%Mass
M1=density.*S1.*Leff1;
M2=density.*S2.*Leff2;
%Matrix element
a11=((-M1*(w.^2))+((density.*c^2.*(S1)^2)./V1));
a12=((-density.*(c^2).*S1.*S2)./V1);
a21=((-density.*(c^2).*S1.*S2)./V1);
a22=(-M2.*(w.^2))+((density.*(c^2)*((S2)^2)*((1./V1)+(1./V2))));
%Results
Z=(((a11.*a22)-(a12.*a21))./(a22.*density.*c.*S1.*w.*1i));
TL=20.*log10(abs(1+(S1./(2.*Sd.*Z))));
%% Optimization function
lb=[10e-5,10e-5,10e-5,10e-5,10e-5,10e-5]; ub=[1,1,1,1,1,1]; x0=(lb+ub)/2;
fun = @(x,f) 20.*log10(abs(1+(x(3)./(2.*Sd.*Z))));
x = lsqcurvefit(fun,x0,f,xn_fft,lb,ub)
plot(f,xn_fft,'ko',f,fun(x,f),'b-')
xlim([0 fNy]); xlabel('Frequency (Hz)'); ylabel('Transmission Loss (dB)');

Accepted Answer

Walter Roberson
Walter Roberson on 14 Jan 2021
clear all;clc;
Okay, so this is a script, not a function.
x=[ ];
Ln1=x(1);Ln2=x(2); % necks length
S1=x(3); S2=x(4); % necks area
V1=x(5); V2=x(6); % Cavities volume
You initialize x to empty, and then try to use 6 values out of it.
  4 Comments
Hasan Hassoun
Hasan Hassoun on 14 Jan 2021
@Walter Roberson'fun' is the file name of the following function
%% The function is:
function TL = six_var(x0,f)
c=340.3; % Speed of sound in the air
density=1.225; % Air density
Sd=4.3*4.3e-4; % duct area
Lv1=2.5e-3;Lv2=2.5e-3; %length of the cicular cavities
w=2*pi()*f;
%Variables
Ln1=x0(1);Ln2=x0(2); % necks length
S1=x0(3); S2=x0(4); % necks area
V1=x0(5); V2=x0(6); % Cavities volume
Rn1= (S1/(pi()))^0.5; Rn2= (S2/(pi()))^0.5; % necks area
Rv1=((V1/(pi()))^0.5)/0.05; Rv2=((V2/(pi()))^0.5)/0.05;
% Calculations
% Effective neck length
correc_n1d=0.82*(Rn1/2);
correc_n1v1=0.85*Rn1*(1-(1.25*(Rn1/Rv1)));
correc_n2v1=0.85*Rn2*(1-(1.25*(Rn2/Rv1)));
correc_n2v2=0.85*Rn2*(1-(1.25*(Rn2/Rv2)));
Leff1=Ln1+correc_n1d+correc_n1v1;
Leff2=Ln2+correc_n2v1+correc_n2v2;
%Mass
M1=density.*S1.*Leff1;
M2=density.*S2.*Leff2;
%Matrix element
a11=((-M1*(w.^2))+((density.*c^2.*(S1)^2)./V1));
a12=((-density.*(c^2).*S1.*S2)./V1);
a21=((-density.*(c^2).*S1.*S2)./V1);
a22=(-M2.*(w.^2))+((density.*(c^2)*((S2)^2)*((1./V1)+(1./V2))));
%Results
Z=(((a11.*a22)-(a12.*a21))./(a22.*density.*c.*S1.*w.*1i));
TL=20.*log10(abs(1+(S1./(2.*Sd.*Z))));
end
Walter Roberson
Walter Roberson on 14 Jan 2021
You can do less-bad by using different x0. But it takes work to get anything useful.
rng(654321)
format long g
% Harmonic sine waves
fs = 400; % Sampling rate [Hz]
Ts = 1/fs; % Sampling period [s]
fNy = fs / 2; % Nyquist frequency [Hz]
duration = 0.5; % Duration [s]
t = 0 : Ts : duration-Ts; % Time vector
noSamples = length(t); % Number of samples
f = 0 : fs/noSamples : fs - fs/noSamples; % Frequency vector
fundamental_signal = 1.0.*sin(2 .* pi .* 50 .* t); %Fundamental signal
First_Harmonic = 1.0.*sin(2 .* pi .* 100 .* t); %FH1=100 Hz %First Harmonic
Total_signal = fundamental_signal + First_Harmonic; % Contaminated signal
xn_fft = abs(fft(Total_signal)/noSamples); % FFT
%% Optimization function
lb=[10e-5,10e-5,10e-5,10e-5,10e-5,10e-5]; ub=[1,1,1,1,1,1]; x0=rand(1,6);
options = optimset('FunValCheck', 'on', 'MaxFunEvals',10e10, 'MaxIter', 1e5, 'Display', 'iter');
[x,resnorm,residual,exitflag,output] = lsqcurvefit(@six_var,x0,f,xn_fft,lb,ub,options)
Norm of First-order Iteration Func-count f(x) step optimality 0 7 342724 9.43e+04 1 14 286019 0.457509 1.34e+05 2 21 237240 0.400185 2.74e+04 3 28 197157 0.542661 4.2e+04 4 35 190392 0.0609414 4.14e+04 5 42 187581 0.0405134 3.46e+04 6 49 181324 0.0606633 7.42e+04 7 56 177100 0.0597402 2.1e+05 8 63 176184 0.015908 3.65e+04 9 70 175652 0.00397701 3.3e+04 10 77 174704 0.00795402 2.7e+04 11 84 172980 0.015908 2.29e+04 12 91 172619 0.0318161 9.47e+05 13 98 169580 0.00795402 3.67e+04 14 105 169068 0.00795402 4.75e+04 15 112 168908 0.00795402 4.89e+04 16 119 168604 0.00198851 2.47e+04 17 126 168238 0.00397701 2.02e+04 18 133 167540 0.00795402 1.7e+04 19 140 166573 0.015908 6.53e+04 20 147 166573 0.0135607 6.53e+04 21 154 166213 0.00339017 1.74e+04 22 161 165918 0.00339017 1.7e+04 23 168 165330 0.00678035 1.74e+04 24 175 164159 0.0135607 1.65e+04 25 182 161858 0.0271214 2.64e+04 26 189 161599 0.0542428 2.22e+06 27 196 157505 0.0135607 1.18e+05 28 203 157267 0.0135607 3.61e+05 29 210 156290 0.00339017 2.11e+04 30 217 156000 0.00339017 1.42e+04 31 224 155428 0.00678035 1.93e+04 32 231 154465 0.0135607 4.28e+04 33 238 152495 0.0271214 6.68e+04 34 245 152495 0.0130586 6.68e+04 35 252 152176 0.00326465 1.53e+04 36 259 151881 0.00326465 1.46e+04 37 266 151336 0.00652931 1.47e+04 38 273 150326 0.0130586 1.87e+04 39 280 149089 0.0261172 1.07e+05 40 287 148302 0.0261172 4.86e+05 41 294 145456 0.00652931 5.07e+04 42 301 143840 0.00652931 4.94e+04 43 308 141039 0.0130586 4.15e+04 44 315 136532 0.0261172 3.39e+04 45 322 129797 0.0522345 3.58e+05 46 329 112305 0.104469 3.85e+04 47 336 85492.6 0.208938 2.51e+05 48 343 72587 0.112338 2.95e+04 49 350 45371.9 0.131663 1.77e+05 50 357 44523.4 0.0175955 9.97e+04 51 364 44523.4 0.00147913 9.97e+04 52 371 44504.4 0.000369783 1.04e+04 53 378 44480.4 0.000369783 9.98e+03 54 385 44436.2 0.000739567 9.18e+03 55 392 44359.6 0.00147913 7.73e+03 56 399 44231.8 0.00295827 6.97e+03 57 406 43991.8 0.00591653 8.37e+03 58 413 43541.6 0.0118331 8.8e+03 59 420 43541.6 0.0236661 8.8e+03 60 427 43458.8 0.00591653 3.19e+05 61 434 43458.8 0.00591653 3.19e+05 62 441 43458.8 0.00147913 3.19e+05 63 448 43342.2 0.000369783 2.62e+04 64 455 43337.7 0.000369783 5.92e+03 65 462 43333.6 9.24458e-05 5.96e+03 66 469 43325.4 0.000184892 6.06e+03 67 476 43309.4 0.000369783 6.26e+03 68 483 43278.4 0.000739567 6.63e+03 69 490 43217.7 0.00147913 7.26e+03 70 497 43097.5 0.00295827 7.97e+03 71 504 42858 0.00591653 7.88e+03 72 511 42415.8 0.0118331 9.01e+03 73 518 42415.8 0.0236661 9.01e+03 74 525 42349.4 0.00591653 3.49e+05 75 532 42349.4 0.00147913 3.49e+05 76 539 42220.8 0.000369783 3.36e+04 77 546 42220 0.000369783 5.64e+03 78 553 42215.7 9.24458e-05 5.68e+03 79 560 42207.4 0.000184892 5.76e+03 80 567 42191 0.000369783 5.93e+03 81 574 42159.8 0.000739567 6.32e+03 82 581 42099.3 0.00147913 7.02e+03 83 588 41979.7 0.00295827 7.86e+03 84 595 41741.8 0.00591653 7.62e+03 85 602 41308.4 0.0118331 9.06e+03 86 609 41107.2 0.0236661 2.64e+05 87 616 38171.9 0.00591653 3.36e+04 88 623 36620 0.00591653 2.09e+04 89 630 34209.1 0.0118331 1.55e+04 90 637 32218.2 0.0236661 1.29e+06 91 644 32153.6 0.00180947 9.47e+03 92 651 32103.3 0.000452368 9.15e+03 93 658 32007.7 0.000904736 8.48e+03 94 665 31838.6 0.00180947 6.84e+03 95 672 31623.4 0.00361894 5.86e+03 96 679 31497.9 0.00361894 5.88e+03 97 686 31282.9 0.00723789 6.26e+03 98 693 31282.9 0.0144758 6.26e+03 99 700 31198.3 0.00361894 1.31e+05 100 707 31198.3 0.00361894 1.31e+05 101 714 31198.3 0.000904736 1.31e+05 102 721 31177.9 0.000226184 5.43e+03 103 728 31170.4 0.000226184 5.45e+03 104 735 31155.6 0.000452368 5.5e+03 105 742 31125.9 0.000904736 5.6e+03 106 749 31066.9 0.00180947 5.8e+03 107 756 30949 0.00361894 6.05e+03 108 763 30713.6 0.00723789 6.11e+03 109 770 30243.5 0.0144758 5.96e+03 110 777 29378.7 0.0289516 2.34e+05 111 784 27979.7 0.0579031 1.05e+04 112 791 22734.4 0.0579031 6.29e+05 113 798 14406.8 0.109425 4.99e+04 114 805 8753.91 0.12306 5.91e+03 115 812 4696.88 0.1854 1.4e+04 116 819 3518.86 0.0509416 8.49e+05 117 826 3518.86 0.000582084 8.49e+05 118 833 3388.62 0.000145521 1.63e+05 119 840 3387.21 0.000145521 1.08e+04 120 847 3386.16 3.63803e-05 8.64e+03 121 854 3384.43 7.27605e-05 5.6e+03 122 861 3381.77 0.000145521 2.54e+03 123 868 3377.53 0.000291042 961 124 875 3370.92 0.000582084 987 125 882 3362.16 0.00116417 5.62e+03 126 889 3349.07 0.00232834 1.04e+03 127 896 3318.71 0.00465667 8.1e+04 128 903 3284.58 0.00931335 702 129 910 3284.58 0.00931335 702 130 917 3257.52 0.00232834 1.16e+04 131 924 3244.89 0.00232834 972 132 931 3220.39 0.00465667 1.25e+05 133 938 3220.39 0.00465667 1.25e+05 134 945 3220.39 0.00116417 1.25e+05 135 952 3220.39 0.000291042 1.25e+05 136 959 3214.82 7.27605e-05 950 137 966 3214.25 7.27605e-05 955 138 973 3213.14 0.000145521 965 139 980 3211.02 0.000291042 981 140 987 3206.91 0.000582084 998 141 994 3198.78 0.00116417 989 142 1001 3183.04 0.00232834 945 143 1008 3155.71 0.00465667 6.62e+03 144 1015 3103.52 0.00931335 802 145 1022 3001.96 0.0186267 7.36e+04 146 1029 2817.1 0.0372534 514 147 1036 2501.27 0.0745068 7.57e+05 148 1043 2210.84 0.0745068 1.18e+03 149 1050 1659.54 0.0745068 1.84e+06 150 1057 1514.75 0.000303399 673 151 1064 1310.24 0.0745068 1.58e+06 152 1071 1226.37 0.000258548 653 153 1078 1224.04 6.4637e-05 642 154 1085 1219.47 0.000129274 620 155 1092 1210.8 0.000258548 576 156 1099 1195.25 0.000517096 487 157 1106 1171.54 0.00103419 308 158 1113 1155.42 0.00206838 7.16e+04 159 1120 1150.97 0.00206838 255 160 1127 1144.67 0.00206838 1.4e+05 161 1134 1144.67 0.00206838 1.4e+05 162 1141 1144.67 0.000517096 1.4e+05 163 1148 1144.67 0.000129274 1.4e+05 164 1155 1141.94 3.23185e-05 307 165 1162 1141.77 3.23185e-05 308 166 1169 1141.44 6.4637e-05 311 167 1176 1140.8 0.000129274 315 168 1183 1139.6 0.000258548 318 169 1190 1137.34 0.000517096 304 170 1197 1133.61 0.00103419 429 171 1204 1126.33 0.00206838 6.63e+03 172 1211 1111.55 0.00413677 237 173 1218 1083.74 0.00827353 3.23e+04 174 1225 1014.7 0.0165471 2.06e+05 175 1232 732.595 0.036127 295 176 1239 548.215 0.0661883 196 177 1246 233.382 0.132377 72.2 178 1253 63.5293 0.103806 39.3 179 1260 22.572 0.0803194 7.15 180 1267 9.79137 0.0413256 3.32e+03 181 1274 7.6226 0.0107216 1.95e+03 182 1281 7.6226 0.0112791 1.95e+03 183 1288 7.48487 0.00281978 391 184 1295 7.32778 0.00563955 189 185 1302 7.00331 0.0112791 9.43e+03 186 1309 6.74463 0.0049156 6.86e+03 187 1316 6.51793 0.00720694 1.61e+04 188 1323 6.38133 0.00330017 9.16e+03 189 1330 6.07504 0.00711624 1.08e+04 190 1337 5.77238 0.00501194 9.66e+03 191 1344 5.19119 0.00851541 8.01e+03 192 1351 4.59365 0.0068332 6.02e+03 193 1358 4.17335 0.004712 2.16e+03 194 1365 4.14928 0.000276379 2.09e+03 195 1372 4.14807 1.39558e-05 2.09e+03 196 1379 4.14807 7.98848e-08 2.09e+03 Local minimum possible. lsqcurvefit stopped because the size of the current step is less than the value of the step size tolerance.
x = 1×6
0.999999999999933 0.999998683230503 0.000100000000024837 0.000367368999991708 0.00010000213295318 0.00469627172794752
resnorm =
4.14806968555364
residual = 1×200
-8.38257435061974e-16 0.0118366301687791 0.0886431735672132 1.23444700673979 1.24169501112485 0.20683585948775 0.0883990937654325 0.0507299592782027 0.033528486970101 0.0240588074699079 0.0182153949395271 0.0143217885742493 0.0115801954676191 0.00956781541476914 0.00804193263823626 0.00685433860546104 0.00590997684886454 0.00514539181580554 0.0045167917939324 0.00399309560138593 0.0035517073727665 0.00317585492055436 0.00285286016698265 0.00257298441854796 0.00232863886175366 -0.497886166743474 0.00192378363733579 0.00175462838792055 0.00160321957173095 0.00146696740958844
exitflag =
2
output = struct with fields:
firstorderopt: 2090.92117485911 iterations: 196 funcCount: 1379 cgiterations: 0 algorithm: 'trust-region-reflective' stepsize: 7.9884762508321e-08 message: '↵Local minimum possible.↵lsqcurvefit stopped because the size of the current step is less than↵the value of the step size tolerance.↵↵<stopping criteria details>↵↵Optimization stopped because the norm of the current step, 7.988476e-08,↵is less than options.StepTolerance = 1.000000e-06.↵↵'
plotit(x, f, xn_fft, fNy)
residuefun = @(x) sum((six_var(x, f) - xn_fft).^2);
[x_from_search, residue_squared] = fminsearch(residuefun, x0)
x_from_search = 1×6
0.348791933503064 0.799471657195371 -1.34230919409332e-06 0.914085410556639 0.811715480969739 1.00111593242472
residue_squared =
0.999969097845243
residue = sqrt(residue_squared)
residue =
0.999984548803252
plotit(x_from_search, f, xn_fft, fNy)
function plotit(x, f, xn_fft, fNy)
%% plotting
c=340.3; % Speed of sound in the air
w=2*pi()*f;
density=1.225; % Air density
Sd=4.3*4.3e-4; % duct area
Lv1=2.5e-3;Lv2=2.5e-3; %length of the cicular cavities
%Variables
Ln1=x(1);Ln2=x(2); % necks length
S1=x(3); S2=x(4); % necks area
V1=x(5); V2=x(6); % Cavities volume
Rn1= (S1/(pi()))^0.5; Rn2= (S2/(pi()))^0.5; % necks area
Rv1=((V1/(pi()))^0.5)/0.05; Rv2=((V2/(pi()))^0.5)/0.05;
% Calculations
% Effective neck length
correc_n1d=0.82*(Rn1/2);
correc_n1v1=0.85*Rn1*(1-(1.25*(Rn1/Rv1)));
correc_n2v1=0.85*Rn2*(1-(1.25*(Rn2/Rv1)));
correc_n2v2=0.85*Rn2*(1-(1.25*(Rn2/Rv2)));
Leff1=Ln1+correc_n1d+correc_n1v1;
Leff2=Ln2+correc_n2v1+correc_n2v2;
%Mass
M1=density.*S1.*Leff1;
M2=density.*S2.*Leff2;
%Matrix element
a11=((-M1*(w.^2))+((density.*c^2.*(S1)^2)./V1));
a12=((-density.*(c^2).*S1.*S2)./V1);
a21=((-density.*(c^2).*S1.*S2)./V1);
a22=(-M2.*(w.^2))+((density.*(c^2)*((S2)^2)*((1./V1)+(1./V2))));
%Results
Z=(((a11.*a22)-(a12.*a21))./(a22.*density.*c.*S1.*w.*1i));
TL=20.*log10(abs(1+(S1./(2.*Sd.*Z))));
plot(f,xn_fft,'k*-',f,TL,'b+-')
xlim([0 fNy]); xlabel('Frequency (Hz)'); ylabel('Transmission Loss (dB)');
legend('Data','Fitted TL')
title('Data and Fitted TL Curve')
end
%% The function is:
function TL = six_var(x0,f)
c=340.3; % Speed of sound in the air
density=1.225; % Air density
Sd=4.3*4.3e-4; % duct area
Lv1=2.5e-3;Lv2=2.5e-3; %length of the cicular cavities
w=2*pi()*f;
%Variables
Ln1=x0(1);Ln2=x0(2); % necks length
S1=x0(3); S2=x0(4); % necks area
V1=x0(5); V2=x0(6); % Cavities volume
Rn1= (S1/(pi()))^0.5; Rn2= (S2/(pi()))^0.5; % necks area
Rv1=((V1/(pi()))^0.5)/0.05; Rv2=((V2/(pi()))^0.5)/0.05;
% Calculations
% Effective neck length
correc_n1d=0.82*(Rn1/2);
correc_n1v1=0.85*Rn1*(1-(1.25*(Rn1/Rv1)));
correc_n2v1=0.85*Rn2*(1-(1.25*(Rn2/Rv1)));
correc_n2v2=0.85*Rn2*(1-(1.25*(Rn2/Rv2)));
Leff1=Ln1+correc_n1d+correc_n1v1;
Leff2=Ln2+correc_n2v1+correc_n2v2;
%Mass
M1=density.*S1.*Leff1;
M2=density.*S2.*Leff2;
%Matrix element
a11=((-M1*(w.^2))+((density.*c^2.*(S1)^2)./V1));
a12=((-density.*(c^2).*S1.*S2)./V1);
a21=((-density.*(c^2).*S1.*S2)./V1);
a22=(-M2.*(w.^2))+((density.*(c^2)*((S2)^2)*((1./V1)+(1./V2))));
%Results
Z=(((a11.*a22)-(a12.*a21))./(a22.*density.*c.*S1.*w.*1i));
TL=20.*log10(abs(1+(S1./(2.*Sd.*Z))));
end

Sign in to comment.

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!