Function inside function using lsqcurvefit
3 次查看(过去 30 天)
显示 更早的评论
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)');
采纳的回答
Walter Roberson
2021-1-14
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 个评论
Walter Roberson
2021-1-14
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
更多回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Least Squares 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!