NI & MatLab miscommunication

[AgateUtane]

Hello!

I am using MatLab for modal test data acquisition and processing. I use a session based interface and data acquisition toolbox; I have connected my laptop to a National Instruments data card (NI-9234: https://www.ni.com/en-gb/support/model.ni-9234.html); I connect my force transducer ( https://www.bksv.com/en/products/transducers/vibration/Vibration-transducers/force-transducers/8230 ) and 3-axial accelerometer ( https://www.pcb.com/products?model=356a01 ) to this data card.

I have found that when I change the sampling frequency, my FRF (frequency response function) shifts along the frequency axis i.e. for different sampling rates my FRF has different response in respect to frequency. It is illustrated in the image below (legend shows different sampling rates). Theoretically, that should not be the case as the same modes must occur at the same natural frequencies (all curves in the image have to overlay). Just to add, sampling rate is the only parameter I change, all other test parameters remain the same.

I have tried several options:

• I checked my script using a Simulink model which generates excitation and response. For simulated data, the FRF does not change depending on the sampling rate, which means my MatLab script and the way I calculate the FRF is correct;
• In the raw test data I noticed inconsistent time step, which means that the actual sampling rate is not what I set it to be and the actual rate is not consistent during the test;
• Even with a re-calculated actual sampling rate I was not able to acquire FRFs that overlap;
• In the NI-9234 data sheet I noticed an equation to calculate the sampling rates that the hardware supports, so I used those in my tests, and that still did not solve the issue.

Has anyone else come across a similar issue, please? Any ideas of how to solve this would be much appreciated!

This is the code I use to acquire and process the test data:

%% %-----SESSION SETUP-----%

clc %clear command window

clearvars %clear all variables from the workspace

close all %delete all figures

% This code is for modal test data acquisition only.

% Direction explanation: X - tangential(T), Y - radial(R), Z - lateral(L).

daq.getDevices %display connected devices

session_1=daq.createSession('ni'); %create session for a specific vendor

fs=input('Choose a sampling rate [51200/25600/12800/10240/6400/5120/3200/2560/2048]: '); %[Hz]

t=60; %test duration [s]

session_1.Rate=fs;

n=fs*t; %total number of data points

%add input channels to the session

addAnalogInputChannel(session_1,'Dev1','ai0','IEPE'); %force transducer

addAnalogInputChannel(session_1,'Dev1','ai1','Accelerometer'); %accelerometer X

addAnalogInputChannel(session_1,'Dev1','ai2','Accelerometer'); %accelerometer Y

addAnalogInputChannel(session_1,'Dev1','ai3','Accelerometer'); %accelerometer Z

%current in Amps used to excite accelerometer/force transducer

session_1.Channels(1).ExcitationCurrent=0.002;

session_1.Channels(2).ExcitationCurrent=0.002;

session_1.Channels(3).ExcitationCurrent=0.002;

session_1.Channels(4).ExcitationCurrent=0.002;

%% %-----VARIABLES-----%

acc_sensitivity=0.005; %accelerometer sensitivity [V/g]

ft_sensitivity=0.1085; %force transducer sensitivity [V/N]

%sets accelerometer sensitivity in V/g

session_1.Channels(2).Sensitivity=acc_sensitivity;

session_1.Channels(3).Sensitivity=acc_sensitivity;

session_1.Channels(4).Sensitivity=acc_sensitivity;

%add output channels to the session

addAnalogOutputChannel(session_1,'Dev3',0,'Voltage');

%% %-----EXCITATION SIGNAL SETUP-----%

pd=makedist('Normal','mu',0,'sigma',1); %pd-probability distribution % create normal distribution with mean of 0 and standard deviation of 1

white_noise= random(pd,[n 1]);

norm_factor=max(white_noise);

output_data=white_noise./norm_factor;

figure(1)

plot(output_data);

title('White noise excitation signal')

xlabel('Data points')

ylabel('Normalised white noise magnitude')

queueOutputData(session_1,output_data); %prepare the output signal for application

%% %-----DATA ACQUISITION-----%

[data,time]=session_1.startForeground(); %start data acquisition

g=9.81; %gravitational acceleration [ms-2]

a=data(:,2:4).*g; %acceleration [ms-2]

f=data(:,1)./ft_sensitivity; %force transducer voltage conversion to force [N]

figure(2)

subplot(2, 2, 1)

plot(time,f);

title('Force transducer signal')

xlabel('Time (s)')

ylabel ('Force (N)')

subplot(2, 2, 2)

plot(time,a(:,1));

title('Accelerometer - X signal')

xlabel('Time (s)')

ylabel ('Acceleration (m/s^2)')

subplot(2, 2, 3)

plot(time,a(:,2));

title('Accelerometer - Y signal')

xlabel('Time (s)')

ylabel ('Acceleration (m/s^2)')

subplot(2, 2, 4)

plot(time,a(:,3));

title('Accelerometer - Z signal')

xlabel('Time (s)')

ylabel ('Acceleration (m/s^2)')

%% %-----TRANSFER FUNCTIONS & COHERENCE-----%

% %actual sampling rate calculation

% time_step=diff(time);

% fs=mean(1./time_step);

%basic accelerance transfer function acquisition (input, output, window, noverlap, nfft, fs)

[acc_X,freq]=tfestimate(f,a(:,1),[],[],[],fs);

[acc_Y,]=tfestimate(f,a(:,2),[],[],[],fs);

[acc_Z,]=tfestimate(f,a(:,3),[],[],[],fs);

direction=input('Choose direction of interest [X/Y/Z]: ','s');

if direction=='X'

acc=acc_X;

elseif direction=='Y'

acc=acc_Y;

else

acc=acc_Z;

end

figure(3)

plot(freq, abs(acc))

str=sprintf('Accelerance for %c direction versus frequency',direction);

title(str)

xlabel('Frequency (Hz)')

ylabel('|Accelerance| (ms^{-2}/N)')

 

[croohcifer]

So, I'm not super familiar with the DAQ toolbox from MATLAB, but I have a guess as to what might be happening.  In LabVIEW/C/C++/C# implementations of the NI-DAQmx drivers, there is a direct call to the driver to query the actual (i.e., coerced) sample rate.  I'm not sure how the "time" vector gets generated (it's don't think it's the NI-DAQmx driver, so I'm guessing it's MATLAB code that does this), but it could be that it's not using that call under the hood.  I'm also not sure exactly how rounding in MATLAB works.  If you, for example, ask the driver to sample at 2048 S/s, but MATLAB represents that value as 25600.00001656845 or whatever, it could be that the DAQmx driver rounds up to the next supported sample rate (51200 S/s).  If MATLAB then decides to generate the time vector based on the assumed sample rate (25600 S/s, let's say) and not the actual sample rate (51200 S/s, let's say), then you could see the issue you're seeing.  The LabVIEW/C/C++/C# APIs from NI take of these niche rounding cases, but I'm not sure if the MATLAB one does.

Basically, I'm not sure how trustworthy the time vector is since it's almost certainly not NI generating it.  (I say that because none of the NI-provided APIs in any language I've used generate a time vector, just a t0 and a dt, which is all the info you need to reconstruct a time vector without wasting RAM to represent every single time stamp for every sample like the DAQ toolbox appears to be doing.)

If you set a sample rate of 2048 and record data for 60 seconds, you should see 2048*60 samples in the returned data and time vector.  If you see more samples, then you know it's probably this rounding error taking place.  If a rounding error is confirmed to be the culprit, one "easy" fix might be to subtract 1 from your selected sample rate and pass that into the DAQ toolbox - I.e., change "session_1.Rate=fs;" to "session_1.Rate=(fs - 1);".  There are no 9234-supported sample rates that differ by only 1 Hz, so this should always be safe for the 9234, and the true sample rate will always be rounded up to the next supported sample rate.  If the time vector is really generated by MATLAB, as I suspect, based on the requested and not the actual sample rate, then it should be unaffected, and the rest of your code's logic should work.

[AgateUtane]

Thank you for your reply!

I checked my data files and I do get the correct number of samples in my returned data and time vectors (sampling rate * duration). There is, however, something weird going on with the time vector as the time step seems to be inconsistent. 

For example, in one of my data files, the time vector looks like this:

 0
4.882800000000000e-04
9.765600000000000e-04
1.464840000000000e-03
1.953120000000000e-03
2.441400000000000e-03
2.929680000000000e-03
3.417960000000000e-03
3.906240000000000e-03
4.394520000000000e-03

The difference between two consecutive time entries (time step) is 

4.882800000000000e-04
4.882800000000000e-04
4.882800000000000e-04
4.882800000000000e-04
4.882799999999998e-04
4.882800000000002e-04
4.882800000000002e-04
4.882799999999998e-04
4.882799999999998e-04
4.882799999999998e-04

Now, theoretically, all time steps should be the same, however, that is not the case. I am not sure what causes MatLab to construct an inconsistent time step, but this is obviously ruining the sampling rate, which in this case should be 2048, but in reality taking 1/time_step, I get 

2.048005242893422e+03
2.048005242893422e+03
2.048005242893422e+03
2.048005242893422e+03
2.048005242893423e+03
2.048005242893421e+03
2.048005242893421e+03
2.048005242893423e+03
2.048005242893423e+03
2.048005242893423e+03

[AgateUtane · Accepted solution]

If anyone was following this post, I investigated this issue further by contacting National Instruments and MathWorks technical support teams. Turns out that there is a bug in MatLab that causes data disruption at certain sampling rates when NI-9234 is used. MathWorks informed me that the development team is investigating this bug now, so hopefully it will get fixed soon!

Thank you!

[Kevin_Price]

FWIW, the so-called differences in the time steps you posted are utterly negligible.  You've got 15 significant digits of total agreement, and it's no surprise at all that tiny discrepancies show up at the 16th -- it's simply the nature of 64-bit floating point representation.

So again, I don't know what's meant about "data disruption" at certain rates.   But I'm still inclined to put more scrutiny on the Matlab code, including some of the things croohcifer already brought up.

I see 3 plots that are very similar in shape but with a different amount of "stretch" along the freq axis.  The farther you move to the right, the more they're shifted relative to one another.  This would be consistent with a small error in the "delta freq" because the visible "shift" increases as the frequency increases.

The 4th plot starts by seeming to have a "delta freq" of about 1/2 what it should be.  But eventually, the calculated response takes on some different shape characteristics than the others.

I've never used the data acq toolbox and have pretty minimal knowledge of Matlab.  But I see you assign an intended acquisition time t and a # samples n, but I don't see them used to define your data acq task.  At present your requested sample rate is input by the operator.  It appears you used to calculate it from an average of the time steps.  I would expect the latter method to be nominally more likely to be accurate.  However, if you didn't clear out your time vector between runs, a later run might not overwrite *all* the time values held over from the prior run.  That could easily lead to enough error in the average time step to explain the "stretch" effect seen in 3 of the plots.

But at the moment, that calculation is commented out.  It just isn't clear whether it was commented out before or after you generated the plots.

Honestly, I'm a little dubious that your issues would be explained by the "data disruption" idea.  That kind of bug would affect your time domain data.  If you lost raw time domain data, I don't think you'd get 3 frequency domain plots that were ready to overlay almost perfectly, subject to a small correction to their respective frequency increments.

Whatever that bug is, I'm still not so sure you can't get reasonable results just by making some adjustments in your Matlab code.

-Kevin P

[AgateUtane]

Hi, Kevin!

Thank you for your comments!

I discussed this issue with MathWorks technical support, and they confirmed that it is indeed a bug, and they are working to solve it. 

Before that, I went through numerous MatLab tutorials about data acquisition, and I modified my code several times trying lots of different settings including setting duration and # of samples as a session property. It made no difference.

MathWorks suggested to test a very simple setup to see if the problem still persists. Basically, I connected my NI-6341 card to my laptop; the card was then directly connected to NI-9234 card, and this card connected back to laptop. I generated a simple sine wave with a frequency of 2 Hz, and I plotted output signal and the input signal read by NI-9234. In this simple test any other test equipment as shakers and transducers were excluded; additionally there were no calculations of transfer functions - I simply plotted two unprocessed signals. 

The code I used is shown here:

clc %clear command window
clearvars %clear all variables from the workspace
close all %delete all figures

daq.getDevices; %display connected devices
session_1=daq.createSession('ni'); %create session for a specific vendor

fs=input('Input a sampling rate in Hz: '); %[Hz]
t=10; %test duration [s]
session_1.Rate=fs;
n=fs*t; %total number of data points

%add input channels to the session
addAnalogInputChannel(session_1,'Dev1','ai0','Voltage'); %excitation signal

%add output channels to the session
addAnalogOutputChannel(session_1,'Dev3','ao1','Voltage');

dt=1/fs;
wave_freq=2;
duration=0:dt:t;
output_data=sin(2*pi*wave_freq*duration);
output_data=output_data';

% figure
% plot(output_data);
% title('Sine wave excitation signal')
% xlabel('Data points')
% ylabel('Sine wave magnitude')

queueOutputData(session_1,output_data); %prepare the output signal for application

[data,time]=session_1.startForeground(); %start data acquisition

measured_signal=data(:,1);

figure
plot(time,output_data,time,measured_signal)
legend

This plot below shows the generated sine wave in blue and the measured signal in orange; sampling rate 3200 Hz. There is a slight lag there but at least the measured signal is also a sine wave of 2 Hz, and it follows the trend of the input signal.

Now, this plot below again shows input signal in blue and output in orange at a sampling rate of 5120 Hz. It is obvious that the measured output signal has some sort of data disruption, which, I believe, causes my transfer function curves to deviate.

I find it very interesting what happens at the sampling rate of 25600 Hz. In my initial post, data measured at this sampling rate deviate the most. Plot below shows that at this sampling rate, the output signal only shows a sine wave with a frequency of ~1 Hz, while the input is 2 Hz. 

Now, work-around for this problem would be to use sampling rates that show no data disruption/deviation. However, I have modified my test setup already to exclude NI-9234, and I am using NI-6341 for analogue output signal and analogue inputs. Using the same MatLab script as in my initial post, there is no data deviation anymore.

Thank you!

[croohcifer]

You literally just proved my comment from above... There is simply some kind of rounding error - probably in the representatation of the double precision float sample rate variable - taking place in the MATLAB wrapper around the DAQ functions.

Let's clarify a few things.

When you refer to a "data disruption", what you're actually referring to is perceived drift.  Your sine wave isn't "disrupted" (this isn't a signal processing term), it was just acquired at the wrong sample rate.  At a sample rate of 3200 S/s, your signals look as though they don't drift.  There appears to be no rounding error in this case - i.e., you ask for acquisition at 3200 S/s, and the DAQ driver actually does that.

At 5120 S/s, they appear to drift, where 10 cycles of the input results in only 9 cycles of the output.  This would be consistent with 5120 getting rounded up to the next supported 9234 sample rate (51200/10 = 5120 and 51200/9 = 5688.88888).  The ratio of those two sample rates is 9:10, which is why you see 9 cycles of output for every 10 of input at this sample rate.

At 26500 S/s, they appear to drift, where 2 cycles of input results in only 1 cycle of output.  This would be consistent with 26500 getting rounded up to the next supported 9234 sample rate (51200/2 = 26500 and 51200/1 = 51200).  The ratio of those two sample rates is 1:2, which is why you see 1 cycles of output for every 2 of input at this sample rate.

As I stated above:  

If you see more samples, then you know it's probably this rounding error taking place.  If a rounding error is confirmed to be the culprit, one "easy" fix might be to subtract 1 from your selected sample rate and pass that into the DAQ toolbox - I.e., change "session_1.Rate=fs;" to "session_1.Rate=(fs - 1);".  There are no 9234-supported sample rates that differ by only 1 Hz, so this should always be safe for the 9234, and the true sample rate will always be rounded up to the next supported sample rate.  If the time vector is really generated by MATLAB, as I suspect, based on the requested and not the actual sample rate, then it should be unaffected, and the rest of your code's logic should work.

Please give that a try and let us know if it solves your problem.

[Kevin_Price]

The analysis from croohcifer sounds spot-on to me.  

See if your Matlab toolbox provides functions for querying the actual sample rate(s) used by the output and input tasks (sessions?).  Querying the actual sample rate would help demonstrate the things croohcifer identified that and also be useful for your further processing, including the possibility that you might need to make your own *correct* time vector.

Doing some further investigation and digging like this seems like a better alternative than being stuck waiting on some kind of bug fix from Matlab.  (While I don't *know*, I wouldn't rule out the possibility that Matlab support could have an inaccurate understanding of NI's products.  What they think is their own bug might be partly due to insufficient knowledge about sample rate quantization on devices like the 9234.)  

-Kevin P

[AgateUtane]

Thank you croohcifer and Kevin for your comments!

I tried setting session sampling rate to (fs-1) however I still see that the output signal does not follow the input signal. The only way I know how to check the sampling rate in MatLab is to calculate 1/time_step and take the mean value of that. I haven't found any other way.

These are the results I got by trying (fs-1):

2048 Hz

5120 Hz

25600 Hz

As I mentioned before, I have excluded NI-9234 from all of my tests and I am using NI USB-6341 instead where I do not get this issue.

Thank you!

[Kevin_Price]

1. It's good to know you've got a workaround using a USB-6341.

2. I tried looking up a function list for Matlab's data acquisition toolbox to see if I could spot something useful.  Wow, the list of available functions is shockingly small!

3. If you've got a trusted function generator and scope, you could make some headway characterizing the issue with the 9234.  That'll be up to you of course, you've got a workaround already.  But the 9234 *does* have some nice features (24-bit conversions, integrated anti-alias filter...) which could make it worthwhile.

    You'd just generate a steady sine or square wave at a known frequency while capturing it with the 9234 at various sample rates.  Then you can post-process the data and times, which might let you figure out the needed compensation factor to cover over the bug in the Matlab toolbox.

-Kevin P