Calculating FRF of signals with disalignment and different sampling frequencies
8 ビュー (過去 30 日間)
古いコメントを表示
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);
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1433583/image.jpeg)
5 件のコメント
Star Strider
2023 年 7 月 16 日
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.
.
回答 (0 件)
参考
カテゴリ
Help Center および File Exchange で Spectral Estimation についてさらに検索
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!