Function inside function using lsqcurvefit
3 views (last 30 days)
Show older comments
Hasan Hassoun
on 14 Jan 2021
Commented: Image Analyst
on 23 Jan 2021
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)');
1 Comment
Accepted Answer
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
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)
plotit(x, f, xn_fft, fNy)
residuefun = @(x) sum((six_var(x, f) - xn_fft).^2);
[x_from_search, residue_squared] = fminsearch(residuefun, x0)
residue = sqrt(residue_squared)
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
More Answers (0)
See Also
Categories
Find more on Least Squares in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!