how to validate mscohere?

I would like to validate the use of mscohere for impact test data, using an instrumented hammer that measures force and an accelerometer to measure response. I ran some artificial cases - see attached pdf. Single input, single output. I tried to force mscohere to use only 3 windows of data, one for each impact test. The first half of results are for sinusoidal input and output, and the second half are for a simulated impact input. There are 2 time domain plots in there. The inputs are the 3 input vectors stacked vertically, and the 3 output vectors stacked vertically. Sample rate is 1000 Hz, with a 4 second window and 0.25 Hz resolution. I am using a rectangular window:

[coherence, frequencies2] = mscohere(coherenceInput, coherenceOutput, rectwin(windowLength), 0, windowLength, sampleRateHz);

The results look about right to me, but I would like to validate it against some known or anlytical results. Are there any suitable examples available?

Also, how do I check that it is getting the absolute basics right - for example, is it using the 3 windows of data that I want it to use?

Answer 1

[William Rose]

@CM,

You say there are two time domain plots, but I see three. Each of the three time domain plots shows an input and an output. I You say "The inputs are the 3 input vectors stacked vertically, and the 3 output vectors stacked vertically." I would not do this, since it creates a signal which may have discontinuities at the boundaries, and those discontinuities in the output would probably look different if the entire concantenated input signal had been run thorugh the system.

To get a coherence estimate with small error bars, I recommend using 8 or more segments, not three. I say that because coherence is estImated from ratios of cross-power and auto-power spectrum estImates. The power spectrum esitmates have variances on the order of (signal power)/sqrt(number of segments). With three segments, the power spectrum estimate variance is on the order of (signal power)/1.73, which can lead to large error bars on the coherence, when using real data.

I am attaching some files of simulated data I made for coherence testing years ago. Each file has two columns: x (input) and y (output). Fs=1000 Hz; number of samples = 1024. The input (x) is the same in all the files. The input is band-limited white noise (30 to 300 Hz). The output equals the input plus uniform noise. The noise is uniform on [-a.+a]. The noise amplitude a is different in ther different files, and the actual random values of the added noise are different. The noise amplituides and files are:

File a=noise half-amplitude

coh1_01xy.txt 0.1

coh1_03xy.txt 0.3

coh1_10xy.txt 1.0

coh1_20xy.txt 2.0

We expect the coherence to be lower when there is more noise on the output, and we exect large error bars on coherence estimates at frequencies where there is little power in the input - because when there is little or no power in the input.

I reocmmend using a Hamming window with half-overlapped segments.

data=load('coh1_10xy.txt'); % produces an error when I run it in Matlab Answers window
                            % but it works when I run it on my computer
x=data(:,1); y=data(:,2);   % x=input, y=output
N=length(data);             % total nmber of samples
Fs=1000;                    % sampling rate (Hz)
Nd=8;                       % number of disjoint segments
nseg=round(N/Nd);           % samples per segment
[cxy,f]=mscohere(x,y,hamming(nseg),round(nseg/2),nseg,Fs);
plot(f,cxy,'-b.'); 
grid on; xlabel('Frequency (Hz)')
title('Coherence'); ylabel('\gamma^2')

Code above produces plot below:

If we plot the coherences for all 4 files, we get

Note that when the noise amplitude (a) is larger, the estimated coherence is lower, as expected. The estimated coherence is low and variable at the frequencies where there is little input power (0 to 30 Hz and >=300 Hz), as expected.

[William Rose]

@CM,

Let X(f) be the input at frequency f, and let Y(f) be the output. Let N(f) be the amplitude of an uncorrelated noise signal, at frequency f.

Suppose , where b and a are real scalars or the magnitudes of complex constants. Then, using the formulas for coherence and variance, one can show that the coherence is

where VarN=|N|^2 and VarX=|X|^2. In the previous example, VarX~=2 between 50 and 300 Hz, VarN~=1/3, b=1, and a=0.1, 0.3, 1.0, or 2.0, for the four files. Therefore the coherences are approximately ~=[0.9983, 0.983, 0.86, 0.60]. The values in the plot above, in the band 50-300 Hz, are around the expected values.

You can use this formula to confirm that mscohere() is doing what you expect, if you make simulated data with known properties.

Remember that we do not measure the true coherence - we just estimate it from our experimental samples. As the plots shows, the coherence estimates are somewhat noisy. They also have bias error. The variance and bias error of the estimate gets smaller when Nd, the number of segments used for estimation, is large.

[CM]

Thanks William. Unfortunately in my case I am geting coherence from impact test data - striking an object with a hammer and measuring the applied force and the dynamic response. So the input and output data cannot be treated as continuous. There has to a discrete number of data blocks analysed, and each block must contain one impact and the 'entire' dynamic response (or at least, until it dies down to the noise floor). So I need to 'trick' mscohere into breaking up the data blocks how I want it to.

Maybe I should be calculating the coherence the long way and comparing that to the mscohere output. That would at least validate that mscohere is breaking up the data as I intended it to. Then I could use your equation to validate some test cases against the theoretical coherence.

Is there an equation for the expected coherence estimate in situations like the one I tested where there is a low number (n) of data windows and the response phase changes by 360/n with each window? Should that give the same coherence as random noise?

Regarding the time waveform plots, I meant there are two plots, each with 3 subplots for the 3 windows of data. Just in case people were confused about what I was using for the input and output.

[William Rose]

@CM,

You can compute coherence from the discrete blows at itrregular intervals. Pick a segment duration that is long enough to have the full impulse and response, with enough time for the system to get back to zero in each segment. It is fine if the exact timing of the impulse is different in different segments. The segments must all have the same duration. THen you may concatenate them all together, as you originally proposed. But if you d this, make sure there is no discontinuity (other than measurement noise) at the segment boundaries. (There shouldn't be discontinuities, because the signals should be zero a the start and end of each segment.)

Since the signal is zero intially and at the end, you don't need a window, i.e. you can use a rectangulkar window which means no window.

I expect the coherence estimates will be quite high, if the different hammer strikes are of similar magnitude, because in that case, I expect the measured responses will be quite similar, and I bexpect measurement noise is not too large. If the strikes vary a lot in magnitude, you might see nonlinearities in the response, i.e. the largre impulses may excite vibration modes that are not seen with weak impulses. Such nonlinear behavior will lead to lower coherence.

Are you doing an experiment similar to what is described here?

If you want to verify that mscohere() is doing the right thing, you could compute the coherence spectrum estimate yourself. But I wouldn;t bother. I did it myself, as a grad student, in BM (before Matlab) era. I don't think it's worth your effort. You'd really just be testing your own code-writing ability against mscohere(), which I would trust more.

But I do recommend that you write code to generate simulated responses to impulses that occur at somewhat random times within the first part of each segment. The simulation should allow for the addition of various amounts of measurement noise. The you can analyze the simulated data with mscohere() and confirm that the results look reasonable. For example, the coherence should be lower when the measurement noise is larger, and when there are nonlinear responses. If the input signal has very little power at some frequencies, then you should see low coherence at those frequencies. The coherence should be unity if there's no noise and if there is some power in the input signal at that frequency.

[CM]

Thanks William. Yes, that example is similar to what I am doing.

I found a bug in my code. I am now consistently getting a coherence of 0 when the 3 responses are offset by 120 degrees at each frequency, a coherence of 0.11111 when the third response is offset by 180 degrees at each frequency, a coherence of 1 (at the 3 tones) when the impact and response are offset by the same amount of time - so the impact occurs at a different location within the window, but still lines up with the same spot on the response function, and a coherence of 0.857 when I swap around the response tone amplitudes of 1, 2 and 3.

[William Rose]

@CM,

Here is a script that simulates repsonse to a sequence of Nd hammer strikes at semi-random intervals. "semirandom" means each strike time is random within the first 0.5 s of each 1 second long window. The strike force is a Gaussian with st.dev.=1 ms. The hammer strike amplitudes are random in the range 100 to 200 N.

The impulse response (acceleration) of the struck object is a set of 8 underdamped resonant oscillations with frequencies approximately 100 Hz to 450 Hz.

The script computes the input force vs. time vector. It computes the output acceleration vs. time in the standard way: by convolving the input with the impulse response. It adds a specified amount of Gaussian random noise to the input and the output. It makes several plots.

hammerImpulseResponse

Resonant frequencies (Hz)=

105.5487 152.7560 199.6176 241.0172 297.9246 355.4221 395.1561 443.2057

After the script finishes, you can save the results to a file:

data=[x;y];
save('hammerData00_00.txt','data',"-ascii")

I used "00_00" in the file name to indicate that the noise amplitude is 0 for the input and 0 for the output, for this example.

You can analyze the simulated data with your coherence script. You can change ANin and ANout in the script, to add different amounts of noise to the input and output, then see how the coherence changes.

[CM]

Thank you very much William. That's a very realistic example.

[CM]

Also, I was able to reproduce the coherence predicted by your equation if, for the incoherent part of the signal, I used 3 sin waves, each offset by 120 degrees. I also tried with random noise, but the results were all over the place, both positive and negative errors of up to 0.7 for the coherence.

I also validated mscohere against the result obtained the 'long' way from cpsd and pwelch, which was the same. I don't think I'll bother looking for the lower lever calculation.

[William Rose]

@CM, good job on calculating the (estimated) coherence spectrum with cpsd and pwelch!