twmacro/pyyeti

IRF through NExT

Closed this issue · 6 comments

Dear @twmacro

I'm curious about how to obtain impulse response functions (IRFs) using the Natural Excitation Technique (NEXT), as the resulting output is a 3D array. Could you help me out here?

Best regards,
Ayubirad

Hi @Ayubirad,

Absolutely! I also plan to add this in the documentation. Using the example that's in the docs up to where it calls ERA, you can call NExT directly:

irf = era.NExT(sol.a, sr, lag_stop=75)
t2 = t[1:76]
for ref in range(irf.shape[2]):
    for out in range(irf.shape[0]):
        plt.plot(t2, irf[out, :, ref])

Does that clear it up?

Best regards,
Tim

Dear @twmacro

Thanks, you see I'm searching for a 2D array that contains the free-decay responses acquired by applying NExT to the ambient vibration data, like IRF of channel 1, channel 2, and so on. Can you please let me know if it's possible here?

Sincerely,
Ayubirad

Hi @Ayubirad,

The 2D IRF array with channel i as the reference is simply irf[:, :, i]. If you want a 2D array of just the IRFs from autocorrelations (assuming all 3 channels were used as reference), you could do this:

i = np.r_[:3]
irf_auto = irf[i, :, i]

Does that clear it up?

Best regards,
Tim

Dear @twmacro

Thank you very much, it worked perfectly! I apologize for getting confused over such a simple matter.
You see, the following code from mathworks fileexchange had me going nuts!:

https://www.mathworks.com/matlabcentral/fileexchange/69510-natural-excitation-technique-next?s_tid=srchtitle

function IRF= NExTT(data,refch,maxlags)

nch=size(data,1);
numref=length(refch);

for refl=1:1:numref            %Loop for reference channels
for chan=1:1:nch               %Loop for all channels

Sensor1=refch(refl);            %No. First sensor to get data from
Sensor2=chan;                  %No. Second sensor to get data from

%--------------------------------------------------------------------------
%Generation of vectors before applying cross-correlation to them

x = data(Sensor1,:);           %Extract data from sensor 1 
y = data(Sensor2,:);           %Extract data from sensor 2

%--------------------------------------------------------------------------
%Generation of cross-correlation vector

XCR=xcorr(x,y,maxlags,'unbiased'); %Cross-correlation Vector with 2*maxlags+1 elements
XCF=XCR(maxlags+1:2*maxlags+1);    %Cross-correlation Vector with maxlags+1 elements

%--------------------------------------------------------------------------
% Generation of cross correlation matrix between reference channels and all
%            channels.its dimensions (nch,numref*n)
           
Y(chan,[refl:numref:refl+numref*(maxlags)])=XCF; %Array for cross-correlation
%--------------------------------------------------------------------------
end
end

IRF=Y;
end

Best regards,
Ayubirad

Hi @Ayubirad,

Copy that! I'm glad you got it working. :-) From a quick look at that code, I definitely prefer the 3D organization of my implementation. I also wonder, is "unbiased" proper? I had thought about that for a while during implementation, but concluded it was correct to NOT use unbiased correlations. I'll go back to the NExT docs to see if they say anything about using the unbiased cross correlations. In any case, in my thinking (which could be wrong), I bet it would only make a very tiny difference to the damping estimates. Still, I'd like to do the correct thing even if the impact is small.

Thanks for the update!

Best regards,
Tim

Hi @Ayubirad,

Just a quick follow-up note. After reviewing the papers and the math more carefully, I'm convinced that using "unbiased" correlations is correct. I've updated NExT in v1.2.9 accordingly and made some other updates as well.

I'll now close this issue ticket.

Thank you,
Tim