for offset = 1:Ls
    r(offset) = 0;
    index1 = offset;
    index2 = offset+Lu;  
    for i = 1: Ns
        r(offset) = r(offset) + r(index1: index1+Lg-1) * r(index2: index2+Lg-1)';
        index1 = index1+Ls;  index2=index2+Ls;  
    end
end

Пусть w=e^(j*2*pi/N)

Например возьмём N=4.
    _               _
   | w^0 w^0 w^0 w^0 |
W= | w^0 w^1 w^2 w^3 |
   | w^0 w^2 w^4 w^6 |
   |_w^0 w^3 w^6 w^9_|
  
X=[x3 x2 x1 x0]' ('- транспонирование)

Обратное дискретное преобразование фурье от X
(масштабирующий коэффициент не рассматриваем):

iDFT(X)=W*X

Последовательности символов в каждой поднесущей:

a0,a1,a2...
b0,b1,b2...
c0,c1,c2...
d0,d1,d2...

Импульсная характеристика сглаживающего фильтра,
например пусть 5 отсчётов на символ, длительность в 3 символа:

h0,h1,h2,h3,h4,h5,h6,h7,h8,h9,h10,h11,h12,h13,h14

Полифазная декомпозиция фильтра:

   _         _
  | h0 h5 h10 |
  | h1 h6 h11 |
H=| h2 h7 h12 |
  | h3 h8 h13 |
  |_h4 h9 h14_|
  
A,B,C,D - векторы столбцы выходы линий задержек с отводами сглаживающих фильтров.

Например:

А0=[0,0,0]'
A1=[a0,0,0]'
A2=[a1,a0,0]'
A3=[a2,a1,a0]'
...

Векторы столбцы выходы чётырёх полифазных фильтров:

H*A,H*B,H*C,H*D

Мультиплексируя в 5 раз чаще символьной частоты с выхода этих
векторов столбцов, от 1-го ряда к 5-му, можем получить четыре
последовательности сглаженных символов,
которые затем нужно подать на комплексные гетеродины.

Выход, вектор столбец, сумма 4-х каналов(.* - поэлементное умножение):
            _   _         _   _         _   _         _   _
OUT1=H*A1.*| w^0 |+H*B1.*| w^0 |+H*С1.*| w^0 |+H*D1.*| w^0 |
           | w^0 |       | w^1 |       | w^2 |       | w^3 |
           | w^0 |       | w^2 |       | w^4 |       | w^6 |
           | w^0 |       | w^3 |       | w^6 |       | w^9 |
           |_w^0_|       |_w^0_|       |_w^0_|       |_w^0_|
            _   _         _   _         _   _         _   _
OUT2=H*A2.*| w^0 |+H*B2.*| w^1 |+H*С2.*| w^2 |+H*D2.*| w^3 |
           | w^0 |       | w^2 |       | w^4 |       | w^6 |
           | w^0 |       | w^3 |       | w^6 |       | w^9 |
           | w^0 |       | w^0 |       | w^0 |       | w^0 |
           |_w^0_|       |_w^1_|       |_w^2_|       |_w^3_|
            _   _         _   _         _   _         _   _
OUT3=H*A3.*| w^0 |+H*B3.*| w^2 |+H*С3.*| w^4 |+H*D3.*| w^6 |
           | w^0 |       | w^3 |       | w^6 |       | w^9 |
           | w^0 |       | w^0 |       | w^0 |       | w^0 |
           | w^0 |       | w^1 |       | w^2 |       | w^3 |
           |_w^0_|       |_w^2_|       |_w^4_|       |_w^6_|
...

Мультиплексируя выход вектора столбца OUT в 5 раз чаще символьной частоты
получим модуляцию FMT в "лоб".

Раскроем OUT1:
      _                                           _
OUT1=| h0*(a0*w^0+b1*w^0+c0*w^0+d0*w^0)+h5*0+h10*0 |
     | h1*(a0*w^0+b1*w^1+c0*w^2+d0*w^3)+h6*0+h11*0 |
     | h2*(a0*w^0+b1*w^2+c0*w^4+d0*w^6)+h7*0+h12*0 |
     | h3*(a0*w^0+b1*w^3+c0*w^6+d0*w^9)+h8*0+h13*0 |
     |_h4*(a0*w^0+b1*w^0+c0*w^0+d0*w^0)+h9*0+h14*0_|

Обозначим F=iDFT([a,b,c,d]')

F1=iDFT([a0,b0,c0,d0]')
F2=iDFT([a1,b1,c1,d1]')
F3=iDFT([a2,b2,c2,d2]')
...

      _                   _
OUT1=| h0*F1(1)+h5*0+h10*0 |
     | h1*F1(2)+h6*0+h11*0 |
     | h2*F1(3)+h7*0+h12*0 |
     | h3*F1(4)+h8*0+h13*0 |
     |_h4*F1(1)+h9*0+h14*0_|
      _                       _
OUT2=| h0*F2(2)+h5*F1(2)+h10*0 |
     | h1*F2(3)+h6*F1(3)+h11*0 |
     | h2*F2(4)+h7*F1(4)+h12*0 |
     | h3*F2(1)+h8*F1(1)+h13*0 |
     |_h4*F2(2)+h9*F1(2)+h14*0_|
      _                           _
OUT3=| h0*F3(3)+h5*F2(3)+h10*F1(3) |
     | h1*F3(4)+h6*F2(4)+h11*F1(4) |
     | h2*F3(1)+h7*F2(1)+h12*F1(1) |
     | h3*F3(2)+h8*F2(2)+h13*F1(2) |
     |_h4*F3(3)+h9*F2(3)+h14*F1(3)_|
...

FF1=[F1 0 0]
FF2=[F2 F1 0]
FF3=[F3 F2 F1]
...

FF - линия задержки с отводами векторами столбцами, матрица 4x3.

OUT1=sum(H.*[FF1(1,:); FF1(2,:); FF1(3,:); FF1(4,:); FF1(1,:)])'
([FF1(1,:); FF1(2,:); FF1(3,:); FF1(4,:); FF1(1,:)] - склейка рядов в матрицу 5x3;
sum() - сумма элементов каждого ряда матрицы, в итоге получаем вектор строку из матрицы)

OUT2=sum(H.*[FF2(2,:); FF2(3,:); FF2(4,:); FF2(1,:); FF2(2,:)])'
OUT3=sum(H.*[FF3(3,:); FF3(4,:); FF3(1,:); FF3(2,:); FF3(3,:)])'
OUT4=sum(H.*[FF4(4,:); FF4(1,:); FF4(2,:); FF4(3,:); FF4(4,:)])'
...

В результате получаем быстрый алгоритм где вместо комплексного гетеродина
на каждую поднесущую имеем быстрое блочное FFT на все поднесущие,
вычисляемое на символьной частоте, вместо полифазного фильтра на каждую
поднесущую имеем один полифазный фильтр на все поднесущие, работающий на
символьной частоте.