Function inside function using lsqcurvefit
2 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Hasan Hassoun
le 14 Jan 2021
Commenté : Image Analyst
le 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 commentaire
Réponse acceptée
Walter Roberson
le 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 commentaires
Walter Roberson
le 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
Plus de réponses (0)
Voir également
Catégories
En savoir plus sur Least Squares dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!