Error in interp1 (line 188). VqLite = matlab.int​ernal.math​.interp1(X​,V,method,​method,Xqc​ol);

function envelope=envfilt(ts,f0,fs)

L=length(ts);

vert=size(ts,2)>size(ts,1);

if vert

ts=ts';

end

inc=find(~isnan(ts));

ts=interp1(inc,ts(inc),(1:L));

W=round(fs/f0);

smoothWindow=normpdf(0:(2*W),W,W/(3*sqrt(2*log(2))))';

odd=mod(W,2);

if odd

e=(1/2)*(imerode(ts,ones(W,1))+imdilate(ts,ones(W,1)));

else

e=(1/2)*(imerode(ts,ones(W,1))+imdilate(ts,ones(W,1)));

e=(1/2)*(e(1:end-1)+e(2:end));

e=interp1((2:L)'-1/2,e,(1:L)','linear','extrap');

end

if vert

envelope=conv(e,smoothWindow,'same')';

else

envelope=conv(e,smoothWindow,'same');

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [tshat,envelope,coeff]=f0param(ts,f0,fs)

%Uses the time series sigma, the guessed frequency f0 and the sample

%frequency fs to parameterize the signal. The output sighat is the

%parameterized reconstruction of sigma, envelope is the envelope of the

%signal (average of sliding minimum and maximum), body is the body

%component of the signal (sliding minimum) and coeff contains the

%coefficients for the sine and cosine functions describing the oscillating

%signal component.

vert=size(ts,2)>size(ts,1);

if vert

ts=ts';

end

L=length(ts);

t=(1:L)'/fs;

envelope=envfilt(sum(ts,2),f0,fs); % Common anvelope for all bands

envelope=envelope/max(envelope); % Same normalization as wave functions

regressor=ones(L,1);

for h=1:floor((fs/2)/f0) % Here it is possible to add harmonics beyound Nyquist

regressor=[regressor sin(2*pi*h*f0*t) cos(2*pi*h*f0*t)];

end

regressor=regressor.*repmat(envelope,[1 size(regressor,2)]); % Apodizing wave functions with envelope (super important step)

for b=1:size(ts,2) %

coeff(:,b)=regressor\ts(:,b);

if vert

tshat(b,:)=(regressor*coeff(:,b))';

else

tshat(:,b)=regressor*coeff(:,b);

end

end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [f0,compcurve,envres,regres,totres,coeff,paramTS,envelope]=find_f0(ts,fs,f0min,f0max,fno,fdetno)

%Finds the frequency that best describes an oscillating signal,

%compensating for biases towards both low and high frequencies. ts is the

%time series, fs is the sample frequency, f0min is the lowest evaluated

%frequency (typically 3*fs/length(ts), which corresponds to the lowest

%frequency at which three full wingbeats fit into ts, but a higher

%frequency may be used), f0max is the highest evaluated frequency

%(typically the Nyquist frequency fs/2, but a lower frequency may be used),

%fno is the number of test frequencies in the coarse evaluation and fdetno

%is the number of frequencies used in the detailed evaluation.

%

%The output f0 is the best frequency the function could find, compcurve is

%the compensated residual curve, envres is the residual of the envelope

%only, regres is the theoretical residual of the regressor based on the

%number of terms included, and totres is the total residual of the time

%series and the reconstruction of the time series.

%

%All residual curves are as function of test frequency, corresponding to

%the vector fTest=linspace(f0min,f0max,fno). The obtained frequency f0

%corresponds to the minimum of the compensated residual compcurve.

%

%The script loops through fTest and evaluates each frequency with the

%function f0param.

fTest=linspace(f0min,f0max,fno); %fTest is vector of frequencies to test

L=length(ts); %L is length of time series

if size(ts,2)>size(ts,1)

ts=ts';

end

%% Loop through test frequencies

for l=1:2

%% Determine mode - coarse or fine

if l==1 %decides if its testing frequency or detailed frequency around frequency

fTestTemp=fTest;

elseif l==2

fTestTemp=linspace(0.9*f0test,1.1*f0test,fdetno);

if fTestTemp(end)>f0max

fTestTemp=fTestTemp(1:find(fTestTemp<f0max,1,'last'));

end

if fTestTemp(1)<f0min

fTestTemp=fTestTemp(find(fTestTemp>f0min,1,'first'):end);

end

end

%% Calibrate frequencies

resEnv=zeros(length(fTestTemp),1); %initiate envelope residual vector

resReg=zeros(length(fTestTemp),1); %initiate regressor residual vector

resTot=zeros(length(fTestTemp),1); %initiate total residual vector

for m=1:length(fTestTemp)

%K=2*floor((fs/2)/fTestTemp(m))+1; %DOF regressor, used to get regressor residual

[paramTS,envelope,coeff]=f0param(ts,fTestTemp(m),fs);

resEnv(m)=sum((ts-envelope*coeff(1)).^2)/L; %residual in each time slot, has to divide by L to get average

resReg(m)=1-length(coeff)/L;

resTot(m)=sum((ts-paramTS).^2)/L;

end

resComp=resTot./(resEnv.*resReg);

[~,f0ind]=min(resComp);

f0test=fTestTemp(f0ind);

if l==1

compcurve=resComp;

envres=resEnv;

regres=resReg;

totres=resTot;

end

end

f0=f0test;

[paramTS,envelope,coeff]=f0param(ts,f0,fs);

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% FUNDAMENTAL FREQUENCIES EXTRACTIONS %%%%%%%%%

clear all

close all

clc

dn='C:\Users\doria\Desktop\CODES\Modification_lidar_code_2023\TAI_MEASUREMENT\';

Fn=dir([dn 'Dualband_Mosquito_CO2_Light_Azh0_StatsMom.mat']);

[~,ind]=sort([Fn.datenum]);

Fn=Fn(ind);

TAB=[];

for k=1:length(Fn)

k

load([dn Fn(k).name]);

TAB=[TAB obs];

end

searchMarg=10;

tp=100;

fs=10000/3;

figure

for k= 1:length(TAB)

k

BS=TAB(k).bsCo';

ts=TAB(k).dt;

L=length(BS);

f0min=3*fs/length(ts);

f0max=fs/2;

F=linspace(fs/L,fs/4,L/2);

W=normpdf(1:(L/1),L/2,L/6);

fbin=linspace(f0min,f0max,round(L/2)+1)';

fbin=(fbin(2:end)+fbin(1:end-1))/2;

body=bodyFit(BS)';

pwrSpec=abs(fft(BS,L));

pwrSpec=2*pwrSpec(1:round(L/2))/L;

[peaks, peakLoc]=findpeaks(pwrSpec);

[pks,~]=max(peaks);

if ~isempty(pks)

f0ind=find(pwrSpec==pks,1,'first');

fno=fbin(f0ind);

fdetno = linspace(fno*(1-searchMarg/100),min(fno*(1+searchMarg/100),f0max),tp);

for i=1:tp

[f0(i),compcurve(i),envres(i),regres(i),totres(i),coeff(i),paramTS(i),envelope(i)]=find_f0(ts,fs,f0min,f0max,fno,fdetno(i)); % THE ERROR IS HERE WHICH CALLS ON THE OTHER THREE FUNCTIONS ABOVE!!!!!!!!!!!

end

[~,ind]=min(compcurve);

TAB(k).f0=f0(ind);

else

close all

TAB(k).f0=NaN;

iter=k;

end

end

For the data it is a 46 MB file whereas the data in Matlabcentral must have a maximum of 5 MB.

[Steven Lord]

What is the full and exact text of the error and/or warning messages you received? Show all the text displayed in orange and/or red in the Command Window. The exact text may be useful in determining what's going on and how to avoid the warning and/or error.

[doria yamoa]

Error using matlab.internal.math.interp1

Interpolation requires at least two sample points for each grid dimension.

Error in inter (line 188)

VqLite = matlab.internal.math.interp1(X,V,method,method,Xqcol);

Error in envfilt (line 14)

ts=inter(inc,ts(inc),(1:L));

Error in f0param (line 22)

envelope=envfilt(sum(ts,2),f0,fs); % Common anvelope for all bands

Error in find_f0 (line 62)

[paramTS,envelope,coeff]=f0param(ts,fTestTemp(m),fs);

Error in exempl_FREQ (line 44)

[f0(i),compcurve(i),envres(i),regres(i),totres(i),coeff(i),paramTS(i),envelope(i)]=find_f0(ts,fs,f0min,f0max,fno,fdetno(i));

Answer 1

[Steven Lord]

So that error is coming from the second of these lines (commented out so that later lines in this answer can run.)

% inc=find(~isnan(ts));
% ts=interp1(inc,ts(inc),(1:L));

What happens if only one value in ts is non-NaN? Then inc will have only one value and in that case the interp1 call will correctly error:

y = interp1(1, 1, 1:2)

Error using matlab.internal.math.interp1
Interpolation requires at least two sample points for each grid dimension.

Error in interp1 (line 188)
VqLite = matlab.internal.math.interp1(X,V,method,method,Xqcol);

since you wouldn't have enough data to do any interpolation. [A different error would occur if none of the values in ts were non-NaN.] I would set a conditional breakpoint on the line where you call interp1 and set the condition that will cause MATLAB to enter debug mode to be numel(inc) < 2. This will let you check that the error only occurs when inc has 1 element.

How to correct the error depends on how ts is created or where it comes from. That's something you'll have to investigate based on your knowledge of the code and the problem the code is trying to solve.

[永新]

出错 interp1 (第 188 行)

VqLite = matlab.internal.math.interp1(X,V,method,method,Xqcol);

出错 wanyou (第 33 行)

Bem=reshape(interp1(tqm,bem,T,'spline'),Y_max*10+1,1);