DWT or MODWT with custom frequency bands
6 件のコメント
回答 (1 件)
Hi @ Michael Hadler,
Addressing your query regarding, “ This lets me guess that the parameters given to 'CustomWaveletFilter' and 'CustomScalingFilter' need to be transformed first to make them useful for the filterbank, but I'm not quite sure how that is achieved. Do you have any ideas?”
By carefully transforming your frequency bands into usable wavelet coefficients and adjusting them according to your sampling rate, you should be able to achieve a more regular response from your custom filter bank. However, I did analyze your code and modified it. The code now displays the default passbands of the filter bank. Passbands are frequency ranges in which the filter bank allows the signal to pass through without significant attenuation.Next, custom frequency bands of interest are defined. These bands represent specific frequency ranges that are of interest for the application at hand. The filter length is determined based on these custom bands, ensuring that it is an even number. Then, it designs a low-pass filter for scaling and a high-pass filter for wavelet using the mean of the custom lower and upper bounds, respectively. The low-pass filter is used to extract the low-frequency components of the signal, while the high-pass filter extracts the high-frequency components. The code checks if the number of elements in each filter is odd. If it is odd, a zero is added to make it even to make sure that the filters have an even number of rows and 1 column, the code checks if the number of elements in each filter is odd. If it is odd, a zero is added to make it even.Finally, a custom filter bank is created using the designed filters. The 'Custom' wavelet option is selected, and the custom high-pass and low-pass filters are specified. The wavelets for both the default and custom filter banks are checked, and the results are plotted for visual verification.The resulting plot shows the time-centered wavelets for both the default and custom filter banks. The wavelets represent the basis functions used in the DWT to decompose the signal into different frequency bands. By visualizing the wavelets, one can gain insights into how the filters are capturing different frequency components of the signal. Here is modified code,
% Define signal properties
signal_length = 4001;
sampling_frequency = 20000;
% Create regular filterbank for reference
max_levels = wmaxlev(signal_length, 'sym4');
fb = dwtfilterbank('SignalLength', signal_length, 'SamplingFrequency', sampling_frequency, 'Level', max_levels);
% Display default bands
bands = dwtpassbands(fb);
disp('Default Bands:');
disp(bands);
% Custom frequency bands of interest
custom_bands = [50 100; 100 300; 300 600; 600 1200];
% Define filter length based on custom bands (ensure it's even)
n = 64; % Ensure this is even
% Low-pass filter for scaling (using mean of custom lower bounds)
lpFilt = fir1(n, (mean(custom_bands(:,1)) / (sampling_frequency / 2)));
% High-pass filter for wavelet (using mean of custom upper bounds)
hpFilt = fir1(n, (mean(custom_bands(:,2)) / (sampling_frequency / 2)), 'high');
% Ensure both filters are column vectors
lpFilt = lpFilt(:); % Ensuring it is a column vector
hpFilt = hpFilt(:); % Ensuring it is a column vector
% Adjust filter sizes to have an even number of rows and 1 column
if mod(numel(lpFilt), 2) == 1
lpFilt = [lpFilt; 0]; % Add a zero to make it even
end
if mod(numel(hpFilt), 2) == 1
hpFilt = [hpFilt; 0]; % Add a zero to make it even
end
% Create custom filter bank using designed filters
fb_custom = dwtfilterbank('SignalLength', signal_length, ... 'SamplingFrequency', sampling_frequency, ... 'Wavelet', 'Custom', ... 'CustomWaveletFilter', hpFilt, ... 'CustomScalingFilter', lpFilt);
% Check wavelets for both filter banks
[psi_default,t_default] = wavelets(fb);
[psi_custom,t_custom] = wavelets(fb_custom);
% Plotting results for visual verification
figure;
subplot(2,1,1);
plot(t_default, psi_default');
title('Default Time-Centered Wavelets');
xlabel('Time');
ylabel('Magnitude');
grid on;
subplot(2,1,2);
plot(t_custom, psi_custom');
title('Custom Time-Centered Wavelets');
xlabel('Time');
ylabel('Magnitude');
grid on;
Please see attached plot.

To address your query regarding,” later cross-correlation analysis (such as modwtcorr)”
Please see updated code below. Add this section to code above to see the plot of the correlation coefficients with confidence intervals.
% Generate two synthetic signals for analysis
t = (0:signal_length-1) / sampling_frequency; % Time vector
signal_1 = sin(2 * pi * 50 * t) + 0.5 * randn(size(t)); % Signal 1: Sinusoidal with noise
signal_2 = sin(2 * pi * 50 * t + pi/4) + 0.5 * randn(size(t)); % Signal 2: Phase-shifted sinusoidal with noise
% Compute MODWT for both signals using the custom filter bank
LEV = max_levels; % Number of levels for MODWT
w_signal_1 = modwt(signal_1, LEV);
w_signal_2 = modwt(signal_2, LEV);
% Display the MODWT coefficients
disp('MODWT Coefficients for Signal 1:'); disp(w_signal_1);
disp('MODWT Coefficients for Signal 2:'); disp(w_signal_2);
% Perform multiscale correlation analysis using modwtcorr function
[wcorr, wcorrci] = modwtcorr(w_signal_1, w_signal_2);
% Display the correlation coefficients and confidence intervals
disp('Correlation Coefficients by Scale:');
disp(wcorr);
disp('Confidence Intervals:');
disp(wcorrci);
%Plot the correlation coefficients with confidence intervals
figure;
plot(wcorr, '-o');
hold on;
errorbar(1:length(wcorr), wcorr, wcorr - wcorrci(:,1)', wcorrci(:,2) - wcorr, 'o');
title('Multiscale Correlation Coefficients with Confidence Intervals');
xlabel('Scale Level');
ylabel('Correlation Coefficient');
grid on;
So, you will see two synthetic signals are generated for analysis, each consisting of a sinusoidal component with noise. The MODWT (Maximal Overlap Discrete Wavelet Transform) is computed for both signals using the custom filter bank. Also, the coefficients for both signals are displayed. The code then performs multiscale correlation analysis using the modwtcorr function, which calculates the correlation coefficients and confidence intervals.Finally, the correlation coefficients with confidence intervals are plotted.Please see attached plot.

4 件のコメント
Hi @Michael Hadler ,
There is no need to sweat. I do understand the frustration. However, please see my response to your recent comments below.
Addressing your query regarding, “Alternatively, when reverting back to the custom filterbank and trying to pass it to modwtcorr as an argument, a similar error occurs (you will have to put the previous section in comments):”
%Modwtcorr with filterbank as argument
[wcorr, wcorrci] = modwtcorr(w_signal_1, w_signal_2,
'FilterBank', fb);
The error message: Error using wavemngr (line 269),Invalid wavelet name: filterbank is indicating that the modwtcorr function is expecting a wavelet name (like 'sym4', 'db1', etc.) rather than a filter bank object. The modwtcorr function does not accept a filter bank directly; instead, it requires a wavelet name to determine the wavelet filters to use. If you refer descriptive section of this function, it will provide you the guidance you seek.
The correct syntax for using modwtcorr according to math works documentation is as follows:
wcorr = modwtcorr(w1, w2, wav);
Where:
*w1 and w2 are the MODWT coefficients of the two signals.
*wav is a string representing the wavelet name.
Here is a complete example that demonstrates how to compute the MODWT correlation between two synthetic signals without encountering the error
% Define signal properties
signal_length = 4001;
sampling_frequency = 20000;
% Create regular filterbank for reference
max_levels = wmaxlev(signal_length, 'sym4');
fb = dwtfilterbank('SignalLength', signal_length,
'SamplingFrequency', sampling_frequency, 'Level', max_levels);
% Generate two synthetic signals for analysis
t = (0:signal_length-1) / sampling_frequency; % Time vector
% Signal 1: Sinusoidal with noise
signal_1 = sin(2 * pi * 50 * t) + 0.5 * randn(size(t));
% Signal 2: Phase-shifted sinusoidal with noise
signal_2 = sin(2 * pi * 50 * t + pi/4) + 0.5 *
randn(size(t));
% Compute MODWT for both signals using the
default wavelet 'sym4'
w_signal_1 = modwt(signal_1, 'sym4');
w_signal_2 = modwt(signal_2, 'sym4');
% Compute wavelet correlation
[wcorr, wcorrci] = modwtcorr(w_signal_1, w_signal_2,
'sym4');
% Display results
disp('Wavelet Correlation Coefficients:');
disp(wcorr);
disp('Confidence Intervals:');
disp(wcorrci);
Going back to, “Is there another way to pass the custom information from the filterbank to these functions?”
Using function handles will allow you to create a wrapper function that will include both the filterbank and any additional parameters you wish to pass. Here’s how you can implement this:
% Define a custom function that takes the filterbank and
additional parameters
function result = customWaveletAnalysis(signal, fb,
customParam)
% Perform MODWT using the provided filterbank
w_signal = modwt(signal, fb);
% Use the custom parameter in some analysis
% For example, scaling the wavelet coefficients
scaledCoefficients = w_signal * customParam;
% Return the result
result = scaledCoefficients;
end
% Example usage
customParam = 2; % Custom parameter to scale the
coefficients
result_signal_1 = customWaveletAnalysis(signal_1, fb,
customParam);
result_signal_2 = customWaveletAnalysis(signal_2, fb,
customParam);
The other option you can use is encapsulating your custom information within a structure which will allow you to pass multiple pieces of related data without altering the function signatures significantly. However, the choice of method will depend on the specific requirements of your analysis and the complexity of the data you wish to pass. By implementing these techniques, you can enhance the flexibility and functionality of your wavelet analysis workflows.
参考
カテゴリ
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!