Calculating FRF of signals with disalignment and different sampling frequencies

Hi all. I've been working on the testing of two different sensors. I would like to calculate the FRF for both sensors using a sine chirp excitation so that I can compare their responses. To compute the FRF I am using an frf.estimator function available in the Mathworks community that uses the pwelch method. However, I can't compute the FRF because of two reasons:

1 - My signals have some disalignment due to the way I acquire them ( it is not possible to do it simultaneously for both sensors and even though I try my best to record at the right timings, I fail everytime, especially because these are sensors with very high sampling rates).

2 - The sensors I need to compare have very different sampling frequencies i.e 26.7kHz and 41kHz. This means that for the same sine chirp excitation, one will get +/-1.65x more samples than the other. I have tried resampling in matlab, however, it seems that it is only possible with integer ratios (2, 3 ,4 etc...).

Below I leave some attempts I made on ways to truncate the signals so that I could align them. I also tried using cross correlation which didn't work very well. Regarding the resampling, I really don't know which approach to take, how can I correctly compare these two sensors? Thank you all very much in advance.

% Data & Sensor Parameters
fs_m = 26700; fs_p = 44100;   % Sampling Frequencies
ts_m = 1/fs_m; ts_p = 1/fs_p; % Sampling Periods

%Truncate response signals based on time domain visualization

truncationDurationBeg = 1; truncationDurationEnd = 1;
truncSampMEMS_beg = round(truncationDurationBeg * fs_m); truncSampPiezo_beg = round(truncationDurationBeg * fs_p); truncSampExc_beg = round(truncationDurationBeg * fs_p);
truncSampMEMS_end = round(truncationDurationEnd * fs_m); truncSampPiezo_end = round(truncationDurationEnd * fs_p); truncSampExc_end = round(truncationDurationEnd * fs_p);
acc_MEMS1 = acc_MEMS1(truncSampMEMS_beg+1:end-truncSampMEMS_end); timeVector1 = timeVector1(truncSampMEMS_beg+1:end-truncSampMEMS_end)-truncationDurationBeg;
acc_Piezo_EU1 = acc_Piezo_EU1(truncSampPiezo_beg+1:end-truncSampPiezo_end); timep= timep(truncSampPiezo_beg+1:end-truncSampPiezo_end)-truncationDurationBeg;
exc = exc(truncSampExc_beg+1:end-truncSampExc_end); t= t(truncSampExc_beg+1:end-truncSampExc_end)-truncationDurationBeg;

% Trying to truncate based on max values ( doesn't seem to work since some
% values in the middle of the signal get eliminated as well)

% thresholdMEMS = max(acc_MEMS1(1:2*fs_m)); % Find value where signal is still "0", this case aprox.1.5 s
% toleranceMEMS = 0.1*thresholdMEMS;  % Tolerance for absolute difference
% transitionIndices_m = abs(acc_MEMS1) > thresholdMEMS + toleranceMEMS; 
% acc_MEMS1 = acc_MEMS1(transitionIndices_m); 
% timeVector1 = timeVector1(transitionIndices_m);
% 
% thresholdPiezo = max(acc_Piezo_EU1(1:2*fs_p)); % Find value where signal is still "0", this case aprox.1.5 s
% tolerancePiezo = 0.1*thresholdPiezo;  % Tolerance for absolute difference
% transitionIndices_p = abs(acc_Piezo_EU1) > thresholdPiezo + tolerancePiezo; 
% acc_Piezo_EU1 = acc_Piezo_EU1(transitionIndices_p); 
% timep = timep(transitionIndices_p);

% Trying Correlation & Alignment

resM = acc_MEMS1; % Response Signal Sensor 1 (Nsamples = 259133)
resP = acc_Piezo_EU1; % Response Signal Sensor 2 (Nsamples = 410130)
resExc = exc; % Excitation Signal (Nsamples = 414540)

% Cross Correlation and Signal Aligment

[C_ExcP, lag_ExcP] = xcorr(resExc, resP); % Longer signal as first argument and the shorter signal as second argument.
[M_ExcP, I_ExcP] = max(abs(C_ExcP));
t_ExcP = lag_ExcP(I_ExcP);

[C_ExcM, lag_ExcM] = xcorr(resP, resM);
[M_ExcM, I_ExcM] = max(abs(C_ExcM));
t_ExcM = lag_ExcM(I_ExcM);

resP = resExc(t_ExcP+1:end);
resM = resExc(t_ExcM+1:end);

[Star Strider]

I am not certain what you want to do. It is certainly possible to use resample using the desired sampling frequency as an input, rather than a ratio. The method I use is:

[sr,tr] = resample(s,t,Fs);

where ‘Fs’ is the desired sampling frequency, ‘s’ and ‘t’ are the original signal and time vector, and ‘sr’ and ‘tr’ are thiir resampled versions. (I would resample both signals to the lower sampling frequency of 26.7 kHz since this avoids the problem of creating data in that signal by upsampling it to the higher sampling frequency.) Name them differently of course, perhaps ‘s1’ and ‘t1’ and ‘sr1’ and ‘tr1’, for one of them, and appropriately different names for the second set.

If you want to model the sensors as linear time-invariant systems, use the resampled signals and the System Identification Toolbox functions (I would start with iddata, then use ssest) and proced from there.

.

[Ivo Bastos]

Thank you @Star Strider. I was a bit confused with the resampling but your way definetely works. My excitation, sensor1 response, and sensor 2 response have, 250980, 259133, 248310 samples, respectively. What I am pursuing is using the pwelch function to calculate PSD and then divide input/output PSD to obtain a FRF for each sensor. I will leave the code below. I was wondering what is you're advice on the best strategy to align the signals in time, and cut some samples out of the larger signals. ( to get same number of samples to allow for the pwelch function to work). Like I said, I've trying using the xcorr function as I showed in the code before but it doesn't produce the desired results. Thank you so much in advance, you have been really helpfull already.

resM = acc_MEMS1; % Response Signal Sensor 1
resP = acc_Piezo_EU1; % Response Signal Sensor 2
resExc = exc; % Excitation Signal

%Resampling exc, MEMS and Piezo
[sm, tm] = resample(acc_MEMS1, timeVector1, fs_m);
[sp, tp] = resample(acc_Piezo_EU1, timep, fs_m);
[sexc, texc] = resample(exc, t, fs_m);

[Sxx,f] = pwelch(res,window,[],nfft,fs); % single-sided PSD
[Sff,~] = pwelch(exc,window,[],nfft,fs); % single-sided PSD
[Sxf,~] = cpsd(res,exc,window,[],nfft,fs); % single-sided CSD ( Cross Power Spectrum)
[Sfx,~] = cpsd(exc,res,window,[],nfft,fs); % single-sided CSD ( Cross Power Spectrum)

switch estimator
    case 'H1'
        H = (Sxf./Sff).'; % H1 estimate ((-) needed??) 
    case 'H2'
        H = (Sxx./Sfx).'; % H2 estimate
    case 'Hv'
        H = ((Sxf./abs(Sxf)).*(sqrt(Sxx./Sff)))';
    otherwise
        error('Unsupported estimator option - in FRFestimate')
end

[Star Strider]

My pleasure!

I am not certain that pwelch is the best option here for the Fourier transforms, with respect to doing the calculations for the FRF determination. The pspectrum function or something like this:

function [FTs1,Fv] = FFT1(s,t)

s = s(:);

t = t(:);

L = numel(t);

Fs = 1/mean(diff(t));

Fn = Fs/2;

NFFT = 2^nextpow2(L);

FTs = fft((s - mean(s)).*hann(L), NFFT)/sum(hann(L));

Fv = linspace(0, 1, NFFT/2+1)*Fn;

Iv = 1:numel(Fv);

FTs1 = FTs(Iv);

end

might be preferable.

Example —

Fs = 100;
T = 15;
t = linspace(0, Fs*T, Fs*T+1)/Fs;
s = sum(sin([1;10;20;30;40]*2*pi*t));

[FTs1,Fv] = FFT1(s, t);

figure
plot(Fv, abs(FTs1)*2)
grid
xlabel('Frequency')
ylabel('Magnitude')

function [FTs1,Fv] = FFT1(s,t)
    s = s(:);
    t = t(:);
    L = numel(t);
    Fs = 1/mean(diff(t));
    Fn = Fs/2;
    NFFT = 2^nextpow2(L);
    FTs = fft((s - mean(s)).*hann(L), NFFT)/sum(hann(L));
    Fv = linspace(0, 1, NFFT/2+1)*Fn;
    Iv = 1:numel(Fv);
    FTs1 = FTs(Iv);
end

Experiment with these to see if you get a more acceptable result.

I will post an Answer when I am satisfied that I have actually provided one. Just now, I am not certain that I have a definitive solution for your problem.

.

[Ivo Bastos]

Thank you once more @Star Strider. Regarding the FRF calculation I am relatively happy with it. The reason for that is that I was given some acquired data both in time domain and frequency domain with the respective FRF calculations that were obtained via a high end software for vibration analysis. Since I was asked to build my own Interface to make these calculations, what I did was sort of "calibrating" my calculations with the ones from the software. I concluded that the pwelch method seems to produce very similar results to the software, which works fine for me. However, the software also did the acquisition of that data, which means that the software ensures the perfect alignment of excitation and response. My biggest challenge has been the resampling and alignment of the signals, as it is an important step to obtain good FRFs. That is why I was asking what would your approach be regarding this issue. At the moment, I am doing it like this, which seems to be closer to what I am looking for. Nevertheless, I still need to find the best way to remove those zeros at the beginning of the signals, and ensure their alignment in a way that would work for different signals ( all of the same type but slighly different start times, number of samples, etc). Thanks in advance!!

%Resampling exc, MEMS and Piezo
[sm, tm] = resample(acc_MEMS1, timeVector1, fs_m);
[sp, tp] = resample(acc_Piezo_EU1, timep, fs_m);
[sexc, texc] = resample(exc, t, fs_m);

[sexc_align, sm_align, Delay] = alignsignals(sexc, sm, 'Method', 'xcorr'); 
sp_align = alignsignals(sp, sm, 'Method', 'xcorr'); 
 
Nsamples_exc = length(sexc_align);
Nsamples_m = length(sm_align);
Nsamples_p = length(sp_align);

diff_mexc = Nsamples_m - Nsamples_exc;
diff_pexc = Nsamples_p - Nsamples_exc;

% Remove samples from the beginning of the signals
sm_align_trimmed = sm_align(diff_mexc+1:end);
sp_align_trimmed = sp_align(diff_pexc+1:end);
sexc_align_trimmed =sexc_align;

% Plot Signals in Time Domain
figure(4);
ax(1) = subplot(3,1,1);
plot(sexc_align_trimmed); xlabel("time(s)"); ylabel("acc (ms-²)"); title("Excitation"); 
ax(2) = subplot(3,1,2);
plot(sp_align_trimmed); xlabel("time(s)"); ylabel("acc (ms-²)"); title("Piezo Signal"); 
ax(3) = subplot(3,1,3);
plot(sm_align_trimmed); xlabel('time(s)'); ylabel('acc (ms-²)'); title('MEMS Signal'); 
linkaxes(ax, 'x');

[Star Strider]

My pleasure!

One way to determine the region where the signals begin is to use the envelope function with the 'peak' argument and a very short window. Then take either the upper or lower envelope results and threshold it to determine where it departs from zero (or some approporiately low value that accounts for any noise in the region up to about 0.5 seconds).

However that may not be necessary if you want the transfer functions (FRFs), since when the begin is less relevant than that the output function slightly lags the excitation function (as it necesarily must). That would appear in the transfer function in the phase respoinse.

I still would not suggest using pwelch to define the FRF vectors here, simply because you actually need the complex Fourier transform vectors to calculate the transfer functions, and pwelch does not provide that information. Using pwelch on the transfer function results would be appropriate (assuming that you need that representation), however calculating the FRF (that I assume are calculated the same way as transfer functions are calculated) needs to be performed first on the complex results. (Ideally the full two-sided fft results could make this more reliable, however using the one-sided results will work.) This means dividing the Fourier transforms of the outputs individually by the Fouirier transform of the input, and plotting that result in the frequency domain, such as I did in my illustration. That will give you the transfer functions (again assuming that it is equivalent to the FRF).

If you then want to invert the transfer functions back to the time domain (use the full, two-sided fft result to do those FRF calculations) and then use pwelch on the inverted result, that would likely work.

.