Need suggestions for calculating the coefficient matrix of Bspline interpolation for any order in MATLAB

21 views (last 30 days)
This algorithm has a complexity of 8K flops (5K multiplications
+ 3K adds) where K is the number of samples,
which is about twice that of the filtering procedure described
in Section 11-C. Moreover, it requires additional
intermediate storage for the modified diagonal elements
of the upper triangular factorization. We note that this
procedure is quite similar to the recursive filtering algorithm
the forward and
backward substitutions correspond to the causal and anticausal
recursive filtering equations, respectively. The only
difference is the use of a nonconstant updating factor ai
in (2.11)a nd (2.12).D uring the forward substitution part
of the algorithm, it is easy to verify that this quantity is given by alpha
Based on this theory need suggestions how could I calculate c(k) alpha_i = -1/(4+alpha(i-1)) initial condition alpha_1= -1/4.This nonlinear difference
equation has a stable steady state Q = -2 + root(3) that
corresponds precisely to the value of the recursive filter
coefficient for the direct cubic spline filter [ 2 ] . In fact, it
can be verified that the convergence of aI to this value is
extremely fast (relative error less than 0.04% error for i
= 3). The same type of analysis can be carried out for the
backsubstitution part of the algorithm. These considerations
clearly show that the filtering and matrix procedures
are essentially equivalent and differ only in their way of
handling the boundary conditions. The former is simpler
and uses fewer operations. Please could you suggest or develop the algorithm based on this mathematical implementation .
I solved the MATLAB code based on the information : Please is it possible to solve :
N = 100;
s = rand(1, N);
% Parameters
iterations = 100; % Number of iterations (or until convergence)
alpha_values = zeros(1, iterations); % Initialize array to store alpha values
% Initial condition
alpha_values(1) = -1/4;
% Iteratively calculate alpha values
for i = 2:iterations
alpha_values(i) = -1 / (4 + alpha_values(i - 1));
end
% Use the last value as the converged alpha
alpha = alpha_values(end);
% Display the result
disp('Alpha values for each iteration:');
disp(alpha_values);
disp(['Converged alpha: ', num2str(alpha)]);
% Value for cubic B-splines
J = min(10, N);
b_i = -alpha / (1 - alpha^2);
% Initialize c_plus array
c_plus = zeros(1, N);
% Step 1: Calculate c^+(1) using Equation (2)
% Ensure we do not exceed the bounds of s by using min(J, N)
c_plus(1) = sum(alpha.^(0:J-1) .* s(1:J));
% Step 2: Calculate c^+(k) for k = 2, ..., N using Equation (3)
for k = 2:N
c_plus(k) = s(k) + alpha * c_plus(k - 1);
end
% Display the result
disp('Initial coefficients c^+(k):');
disp(c_plus);
% Calculate c^-(N) using Equation (4)
c_minus = zeros(1, N);
c_minus(N) = b_i * (2 * c_plus(N) - s(N));
% Calculate c^-(k) for k = N-1, ..., 1 using Equation (5)
for k = N-1:-1:1
c_minus(k) = alpha * (c_minus(k + 1) - c_plus(k));
end
% Display the result
disp('Backward coefficients c^-(k):');
disp(c_minus);
% Calculate final B-spline coefficients c(k) using Equation (6)
c = 6 * c_minus;
% Display the result
disp('Final B-spline coefficients c(k):');
disp(c);
The solution are :
>> SplineCalculation
Alpha values for each iteration:
Columns 1 through 10
-0.2500 -0.2667 -0.2679 -0.2679 -0.2679 -0.2679 -0.2679 -0.2679 -0.2679 -0.2679
Columns 11 through 20
-0.2679 -0.2679 -0.2679 -0.2679 -0.2679 -0.2679 -0.2679 -0.2679 -0.2679 -0.2679
Columns 21 through 30
-0.2679 -0.2679 -0.2679 -0.2679 -0.2679 -0.2679 -0.2679 -0.2679 -0.2679 -0.2679
Columns 31 through 40
-0.2679 -0.2679 -0.2679 -0.2679 -0.2679 -0.2679 -0.2679 -0.2679 -0.2679 -0.2679
Columns 41 through 50
-0.2679 -0.2679 -0.2679 -0.2679 -0.2679 -0.2679 -0.2679 -0.2679 -0.2679 -0.2679
Columns 51 through 60
-0.2679 -0.2679 -0.2679 -0.2679 -0.2679 -0.2679 -0.2679 -0.2679 -0.2679 -0.2679
Columns 61 through 70
-0.2679 -0.2679 -0.2679 -0.2679 -0.2679 -0.2679 -0.2679 -0.2679 -0.2679 -0.2679
Columns 71 through 80
-0.2679 -0.2679 -0.2679 -0.2679 -0.2679 -0.2679 -0.2679 -0.2679 -0.2679 -0.2679
Columns 81 through 90
-0.2679 -0.2679 -0.2679 -0.2679 -0.2679 -0.2679 -0.2679 -0.2679 -0.2679 -0.2679
Columns 91 through 100
-0.2679 -0.2679 -0.2679 -0.2679 -0.2679 -0.2679 -0.2679 -0.2679 -0.2679 -0.2679
Converged alpha: -0.26795
Initial coefficients c^+(k):
Columns 1 through 10
-0.0726 0.7997 0.1755 0.1947 0.3517 0.0022 0.1314 0.9068 0.7131 0.3841
Columns 11 through 20
-0.0431 0.2463 0.2872 0.7443 -0.1840 0.0923 0.1442 0.6105 0.5681 0.4955
Columns 21 through 30
0.3182 0.4618 0.1726 0.6984 0.0018 0.6863 -0.0004 0.3686 0.5269 0.6391
Columns 31 through 40
-0.0901 0.9535 0.5202 0.3474 0.3428 0.3549 0.2112 0.4519 0.3897 0.7132
Columns 41 through 50
0.6037 0.4826 0.2493 0.7448 0.3333 0.2614 0.8690 0.6431 0.3778 0.5212
Columns 51 through 60
0.4474 0.0879 0.2777 0.3965 0.1242 0.8110 -0.0225 0.2320 0.1086 0.1986
Columns 61 through 70
0.3825 0.2086 0.8675 0.1978 0.1318 0.8696 0.7468 0.2388 0.0471 0.2454
Columns 71 through 80
0.3430 0.5030 0.1274 0.5687 0.5588 0.0720 0.0981 0.2704 0.2463 0.3582
Columns 81 through 90
0.4119 -0.0248 0.2691 0.7289 -0.1661 0.9734 0.4695 0.3628 0.4813 0.1083
Columns 91 through 100
0.4298 0.8479 0.3196 0.4355 0.1149 0.4581 0.5013 0.5448 0.2495 0.3006
Backward coefficients c^-(k):
Columns 1 through 10
-0.0740 0.2036 0.0397 0.0273 0.0928 0.0055 -0.0183 0.1997 0.1615 0.1104
Columns 11 through 20
-0.0278 0.0608 0.0195 0.2143 -0.0557 0.0238 0.0035 0.1313 0.1204 0.1190
Columns 21 through 30
0.0516 0.1257 -0.0075 0.2006 -0.0502 0.1893 -0.0202 0.0748 0.0893 0.1938
Columns 31 through 40
-0.0840 0.2235 0.1194 0.0748 0.0684 0.0874 0.0286 0.1044 0.0624 0.1567
Columns 41 through 50
0.1285 0.1241 0.0194 0.1770 0.0843 0.0186 0.1918 0.1530 0.0722 0.1084
Columns 51 through 60
0.1167 0.0117 0.0440 0.1133 -0.0265 0.2230 -0.0211 0.0562 0.0222 0.0256
Columns 61 through 70
0.1031 -0.0023 0.2172 0.0567 -0.0140 0.1840 0.1829 0.0642 -0.0010 0.0507
Columns 71 through 80
0.0561 0.1338 0.0037 0.1134 0.1453 0.0165 0.0104 0.0595 0.0484 0.0657
Columns 81 through 90
0.1130 -0.0099 0.0121 0.2239 -0.1066 0.2317 0.1088 0.0636 0.1254 0.0134
Columns 91 through 100
0.0582 0.2125 0.0547 0.1154 0.0050 0.0963 0.0987 0.1329 0.0488 0.0675
Final B-spline coefficients c(k):
Columns 1 through 10
-0.4441 1.2219 0.2382 0.1638 0.5567 0.0330 -0.1099 1.1983 0.9691 0.6623
Columns 11 through 20
-0.1671 0.3647 0.1171 1.2861 -0.3341 0.1429 0.0208 0.7879 0.7222 0.7137
Columns 21 through 30
0.3093 0.7544 -0.0450 1.2037 -0.3014 1.1358 -0.1209 0.4491 0.5355 1.1625
Columns 31 through 40
-0.5042 1.3411 0.7162 0.4485 0.4105 0.5246 0.1718 0.6261 0.3746 0.9400
Columns 41 through 50
0.7711 0.7446 0.1163 1.0618 0.5058 0.1119 1.1511 0.9178 0.4332 0.6503
Columns 51 through 60
0.7004 0.0705 0.2643 0.6800 -0.1587 1.3378 -0.1266 0.3372 0.1334 0.1535
Columns 61 through 70
0.6186 -0.0139 1.3034 0.3404 -0.0839 1.1040 1.0973 0.3854 -0.0058 0.3045
Columns 71 through 80
0.3363 0.8026 0.0225 0.6807 0.8719 0.0991 0.0621 0.3569 0.2904 0.3941
Columns 81 through 90
0.6781 -0.0595 0.0728 1.3432 -0.6395 1.3900 0.6526 0.3817 0.7522 0.0805
Columns 91 through 100
0.3493 1.2752 0.3284 0.6921 0.0299 0.5778 0.5923 0.7975 0.2927 0.4048
>> Please could you suggest any changes of modification for the code?

Accepted Answer

Shishir Reddy
Shishir Reddy on 6 Nov 2024 at 6:17
Hi Rohitashya
After reviewing your MATLAB code and the algorithm, I have identified a few areas where you can consider optimization:
Vectorization - Utilizing vectorized operations can improve performance by reducing the overhead of loops. So, you can calculate 'c_plus' using vectorized operations.
% Calculate c^+(1)
c_plus(1) = sum(alpha.^(0:J-1) .* s(1:J));
% Calculate c^+(k) for k = 2, ..., N, using vectorized operations
c_plus(2:N) = s(2:N) + alpha * c_plus(1:N-1);
Backward substitution - While direct vectorization of the backward loop is challenging due to its recursive nature, it can still be optimized using MATLAB's 'cumsum' and 'flip' functions to streamline the calculation of 'c_minus'.
c_minus(N) = b_i * (2 * c_plus(N) - s(N));
reversed_c_plus = flip(c_plus(1:N-1)); % Reversing the array
reversed_c_minus = flip(c_minus(N));
cumulative_sum = cumsum(reversed_c_plus);
c_minus(1:N-1) = flip(alpha * cumulative_sum + reversed_c_minus);
c_minus = c_minus - alpha * c_plus;
Although 'cumsum'is typically used for cumulative sums, it can be adapted for some recursive operations by reversing the order of operations.
For more information, regarding the 'cumsum' and, 'flip' functions in MATLAB, kindly refer the following documentations –
I hope this helps.

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!