Discrete Weighted Transform - A Primitive Nth Root of Unity

In summary, this is a conversation about the implementation of the "Discrete Weighted Transform" algorithm, specifically in regards to determining a primitive Nth root of unity in the appropriate domain. The weight signal and weight components are discussed, along with the use of the fast Fourier transform and discrete Fourier transform. Suggestions are made to seek help on a computer science forum and reference is given to a research paper on the topic. A numerical example in section 6 is mentioned as well as steps 2, 3, and 4 in Algorithm W which require a primitive Nth root of unity. A python implementation of the algorithm is also provided.
  • #1
Hello, I got pointed to this forum by OfficeShreder. I have a question I've been puzzling myself over for a while.

I am currently trying to implement the "Discrete Weighted Transform". I have reached a step where I need to determine "a primitive Nth root of unity in the appropriate domain".

I have been reading up on what roots of unity are, but I do not know what the "appropriate domain" is, nor how to determine 'g' according to this.

The research paper can be found here: http://faginfamily.net/barry/Papers/Discrete Weighted Transforms.pdf

The following formula is given:

[tex]X_{k} = \sum_{j=0}^{N-1} a_{j}x_{j}g^{-jk}[/tex]

Along with that formula, in section (2.5), variable 'g', is a primitive Nth root of unity in the appropriate domain.

This is all part of the variant of Algorithm W, as referenced in section 6.
Algorithm W can be found at the end of section 3.

If someone could help me figure out how to determine variable 'g', I would be very appreciative. Thanks in advance.
Last edited:
Mathematics news on Phys.org
  • #2
I'm trying to avoid complaining, because I understand that I have no right to receive your help; it is merely my privilege. But I'm a little stressed out and frustrated as I've posted this on two separate forums and have received very little help (and in this case, none at all). Also, this is fairly time sensitive, and therefore adamently require some assistance.

If you don't know the answer that's fine, but would you be able to point me to another forum or person who may be able to help me?

Sorry if I'm sounding like a jerk in any way. >_>
Last edited:
  • #3
Are you familiar with the usual discrete Fourier transform? This weighted version appears to be a modification the same thing. "appropriate domain" can be whatever ring you are considering your x_i's to be a part of. You could take the complex numbers for example, a primitive nth root of unity would then be e^(2*pi*i/n). You could take the integers modulo something, they talk about number theoretic transform version in section 8 there. What domain you want to take will depend on what you are trying to do. It would probably be fine for you to just work with the complex numbers to get the general idea of what's going on.

You also might have better luck getting help in a computer science forum (sorry don't know any) where people are more likely to know about algorithms like this (I can only claim a passing familiarity). Perhaps the mersenneforum.org people would be a good place to go, they're interested in multiplying large numbers quickly, so will be better acquainted with methods like this.
  • #4
You mean the usual "fast Fourier transform"? Then, no, I have not implemented the FFT before.

My x_i's are determined from the equation in section (6.5):

[tex]x = \sum_{j=0}^{N-1} x_{j}2^{Ceiling(qj/N)}[/tex]

I don't know what "ring" this is part of. But what I am trying to do is what's explained in section 6, according to a variant of Algorithm W from the end of section 3 (as I've stated in the first post).

[tex]n = x^{2} (\bmod{p})[/tex]

I've tried posting this on a computer science algorithm forum, but still have no replies (three forums I've posted this on now). If I keep running out of luck, I'll try the fourth forum you've pointed me to.
Last edited:
  • #5
The fast Fourier transform is a numerical algorithm. shmoe is not talking about that.

Shmoe is talking about the discrete Fourier transform, or DFT, which is essentially the same thing as the normal continuous Fourier transform, with sums in place of integrals.

- Warren
  • #6
Ahh, okay. Well I haven't implemented the DFT either. I'm only working on the DWT.

The reason that confused me is because the DWT is supposedly analogous to the FFT.
  • #7
the FFT is a way to compute the DFT with less operations than the 'naive' way of perform the transformation (just computing the terms in the sum as they appear). the DFT is a special case of the DWT, taking the weights to be 1's.

For their algorithm in section 6, use a primitive Nth complex root of unity. Your weights are irrational, so you can consider your "appropriate domain" to be the complex numbers. They give some details of a numerical example in section 6. Try to fill in the steps with exp^(i*2*Pi/N)=i, as N=4 in this example.
  • #8
Where is the "numerical example in section 6" that demonstrates the step involving finding a primitive Nth root of unity?

I see an example of the entire algorithm at the end, but it skips steps 2, 3 and 4 in Algorithm W. Those steps being the ones that require a primitive Nth root of unity.

It might help slightly to show you where I'm at, with a Pythonic implementation of the algorithm.

import math

A Python Implementation of the IBDWT algorithm, for transferring to C.

    Implementation - Aaron Voelker                  - 28/10/06
    Algorithm      - Richard Crandall & Barry Fagin - 1994

Algorithm W variant: Weighted convolution algorithm for x^2 modulo (2^q-1).

    1) Choose run length N >= q and establish bit-sizes b_j for digits according to (6.6).
       Represent x = {x_j}, according to (6.5), zero-padding only to N digits inclusive.
       Compute the components of the weight signal a according to (6.7).
    2) Compute X = DWT (N, a)x.
    3) Compute Z = X^2.
    4) Compute z = DWT^-1 (N, a)Z.  (This is the weighted convolution x*^a x.)
    5) Set z = Round(z), if noninteger (e.g., floating-point) FFTs were used.
    6) Adjust the digits {z_n} to the digit representation of choice.


    Perform x^2 (mod p).

        x = 78314567209,    q = 37,    p = 2^37-1,    N = 4

        b = {10, 9, 9, 9}
        x = {553, 93, 381, 291}
        a = {1.0, 1.681792830507429, 1.4142135623730951, 1.189207115002721}
        z = {704383., 324600., 523365.0000000001, 463577.9999999999}

        x^2 (mod p) = 704383 + 324600 * 2^10 + 523365 * 2^19 + 463578 * 2^28


def ceiling (z):
    """Round z up to the nearest integer."""
    # See section (3.8).
    # The paper describes the "ceiling" function notation as int (z + 1).
    # However, it's actually supposed to be (z + n), where (n | lim -> 1).
    # Fortunately, the math.ceil function handles this accordingly.
    return int (math.ceil (z))

class algorithm_w:
    def __init__ (self, q, N):
        """Initiate the handle for performing a Discrete Weighted Transform as according to Algorithm W."""
        # See section (6) and end of section (3).
        # This handle will be used to execute Algorithm W, where (x^2 modulo (2^q-1)) is calculated.
        # q denotes that the result will be modulo (2^q - 1).
        # N is the run time of x^2 modulo (2^q - 1).
        self.q    = q
        self.intN = N        # Store as an integer for use as a counter.
        self.N    = float(N) # Store as a float for sufficient accuracy.

    '''Debugging functions.'''

    def sigma (self, x_j):
        """Calculate the value of x based on its base representation, x_j."""
        # See section (6.5).
        # Take the sum of each part times its base.
        # I swear, it's only for debugging purposes, really.
        return sum (x_j[j] * self.x_j_base (j) for j in range (self.intN))

    '''Worker functions.'''

    def qjN (self, j):
        """This little formula is ugly, we'll just push it off to the side."""
        return j * self.q / self.N
    def x_j_base (self, j):
        """Calculate the base of the j'th number in the x series."""
        # See section (6.5).
        return 2 ** ceiling (self.qjN (j))

    '''Algorithm W exclusive functions. End of section (3).'''

    def b_j (self):
        """Calculate the series b_j: The number of bits allocated for digit x_(j-1)."""
        # See section (6.6).
        # Step 1 of Algorithm W.
        # Ugh, why the **** do we need the bit sizes again?
        return [ceiling (self.qjN (j)) - ceiling (self.qjN (j-1)) for j in range (1, self.intN + 1)]

    def a_j (self):
        """Calculate the series a_j: The constant signal, required for a properly weighted convolution."""
        # See section (6.7).
        # Step 1 of Algorithm W.
        return [2 ** (ceiling (self.qjN (j)) - self.qjN (j)) for j in range (self.intN)]

    def x_j (self, x):
        """Calculate the series x_j: The base representation of x, it being the sum of its parts."""
        # Where's the section for this god damnit? Why did I have to figure it out myself?
        # Step 1 of Algorithm W.
        # Divide and modulate in the same manner of how we'd jump from one number system. to another
        # Do this by calculating the base for each part, effectively reversing the sigma.

        # ****, this is probably pretty slow for large values of x. Is there a faster, synthetic way of determining the base representation of x?
        # There has got to be, or else this defeats the whole purpose.

        x_j = []
        for j in range (self.intN):       # Start at the end, since we are reversing the sum.
            base = self.x_j_base (self.N - j - 1)
            x_j  = [int (x / base)] + x_j # Convert from long to integer, and add to the result list.
            x   %= base                   # Determine the newly diminished sum of the sigma notation.
        return x_j

    def dwt (self, x_j, a_j):
        """Calculate X_k = DWT(N, a_j)x_j: The Discrete Weighted Transform, X_k, of x_j according to weight signal a_j, comprised of N non-zero constants."""
        # See section (2.5).
        # Step 2 of Algorithm W.


def run_test ():
    x1   = 78314567209          # this number is supposedly represented as
    x_j1 = [553, 93, 381, 291]  # these four values

    q = 37                      # when these values are true
    N = 4                       # how the hell did we get this number?

    dwt_h = algorithm_w (q, N)

    bits = dwt_h.b_j ()         # calculate the bit sizes for x_j
    print bits

    signals = dwt_h.a_j ()      # calculate the constant signal a
    print signals

    x2 = dwt_h.sigma (x_j1)     # x_j1 produces the x value: 78314567209
    print x1
    print x2
    print x1 == x2              # True

    x_j2 = dwt_h.x_j (x1)       # x1 produces the x_j values: [553, 93, 381, 291]
    print x_j1
    print x_j2
    print x_j1 == x_j2          # True

if __name__ == '__main__':
Last edited:
  • #9
Sane said:
Where is the "numerical example in section 6" that demonstrates the step involving finding a primitive Nth root of unity?

I see an example of the entire algorithm at the end, but it skips steps 2, 3 and 4 in Algorithm W. Those steps being the ones that require a primitive Nth root of unity.

That's the example I was talking about, I didn't say they showed you what primitive root to use or how to find it. I said you can fill in their steps using the primitive root I suggested, i.e. take the data from (6.10) and you can reproduce all the numbers in their example when squaring 78314567209 (mod p) by taking exp^(i*2*Pi/4)=i to be your primitive 4th root of unity. It's been a couple of weeks but I recall it worked out (can't seem to find the scrap paper anymore).

All you were missing was the primitive root, no? Or are there other parts of the algorithm you don't follow? (I haven't read your python code)

1. What is a Discrete Weighted Transform?

A Discrete Weighted Transform (DWT) is a mathematical operation used in signal processing to transform a discrete signal into a different representation. It involves multiplying the signal by a set of weights, which are typically chosen to emphasize certain features or frequencies of the signal.

2. What is a Primitive Nth Root of Unity?

A Primitive Nth Root of Unity is a complex number that, when raised to the power of N, equals 1. It is called "primitive" because it cannot be written as a power of any other root of unity. It is often used in DWT to determine the weights used in the transformation.

3. How is DWT related to Fourier Transform?

DWT is closely related to Fourier Transform, as both are used to analyze signals in terms of their frequency components. However, DWT operates on a discrete signal, while Fourier Transform works on a continuous signal. DWT also uses a different set of basis functions compared to Fourier Transform.

4. What are the applications of DWT?

DWT has many applications in signal processing, including image and sound compression, denoising, and feature extraction. It is also used in data analysis and pattern recognition.

5. How do I calculate the weights for DWT using Primitive Nth Roots of Unity?

The weights for DWT can be calculated using the formula wn = e-2πik/n, where n is the size of the signal and k is an integer between 0 and n-1. This formula uses the primitive nth root of unity, which can be found by solving the equation zn = 1 for z. These weights can then be used in the DWT formula to transform the signal.

Suggested for: Discrete Weighted Transform - A Primitive Nth Root of Unity