Results 1 to 2 of 2

Thread: Here's my own frequency domain transform.

  1. #1

    Thread Starter
    Frenzied Member
    Join Date
    Oct 2008
    Posts
    1,181

    Here's my own frequency domain transform.

    I believe that this transform may be completely new, an invention of my own. It is a variant of the DCT, that I've never seen documented before. It solves a problem that normal DCTs have. It allows you to actually work up to the highest frequency possible. With a normal DCT (lets use an 8-point DCT for this example), when you perform the forward DCT, you have the highest frequency bin representing 7.5 half-cycles per 8-point data set. Yet in an 8 point data set, the highest possible frequency is in fact 8 half-cycles. So while this transform is capable of capturing enough data to be completely reversible and get back your original time-domain data if your starting data is time-domain data, if you are starting in the frequency-domain, and plan to use an IDCT to construct a time-domain signal, you will run into a problem. That problem is that if you set the highest bin to some value in your frequency-domain, and then perform an IDCT, you won't get a nice continuous-amplitude waveform at the highest frequency. Instead you will get a waveform that starts at zero amplitude, and then swells to maximum amplitude in the middle of the data set, and then decreases to zero at the other end. It will be a signal that has the highest frequency, but it's certainly does not have constant amplitude.

    Now that is not fixed in any of the other DCTs, (there's 3 variants, DCT1 which is its own inverse, DCT4 which is its own inverse,and DCT2 which has DCT3 as its inverse). To get around this problem, you can use a DFT where the the bin at datasetsize/2 represents the highest frequency, but for real-valued input the second half of data set is redundant. In the the real-number bins, the second side is a mirror image of the first side. In the imaginary-number bins, the second side is a negative mirror image of the first side.

    This is why I designed the DCT_B (which stands for Discrete Cosine Transform, Ben's version), and its inverse, the IDCT_B. The frequency-domain dataset size is only one bin more than the time-domain dataset size (for any TDSize, FDSize=TDSize+1). So an 8-point time-domain dataset has a 9-point frequency-domain dataset. And the practical effect of this is that when you set the highest frequency bin in the frequency-domain, and then perform an IDCT_B on it, you will get a constant-amplitude signal at the highest frequency.

    Also, whereas other DCTs don't usually include a scaling factor in their definition, and instead choose to describe it as an additional step to normalize the values, the definition of my own DCT, the DCT_B, has a built-in scaling factor. This scaling factor guarenties that if you have a cosine wave, the amplitude of that wave can be read-out directly in the DCT. For example, if I have an 8-point data set containing cos(Pi/8*F*n)*A , where n is a value incremented from 0 to 7 in a for-next loop, and A is the amplitude, and F is an integer between 0 and 8; and then I perform a DCT_B on that data set (which will give me frequency-domain data with bins from 0 to 8), you will find that bin# F contains the value A.

    Here is my VB6 code that defines the DCT_B and its inverse, the IDCT_B.
    Code:
    Private Const Pi As Double = 3.14159265358979
    
    Public Sub DCT_B(ByRef Src() As Double, ByRef Dest() As Double)
        Dim n As Long
        Dim k As Long
        Dim TDLen As Long
        Dim FDLen As Long
        
        TDLen = UBound(Src) + 1
        FDLen = TDLen + 1
        ReDim Dest(FDLen - 1)
        
        For k = 0 To FDLen - 1
            For n = 0 To TDLen - 1
                Dest(k) = Dest(k) + Src(n) * Cos(Pi / TDLen * n * k) / TDLen * 2
            Next n
        Next k
        Dest(0) = Dest(0) / 2
        Dest(FDLen - 1) = Dest(FDLen - 1) / 2
    End Sub
    
    
    
    Public Sub IDCT_B(ByRef Src() As Double, ByRef Dest() As Double)
        Dim n As Long
        Dim k As Long
        Dim TDLen As Long
        Dim FDLen As Long
        
        FDLen = UBound(Src) + 1
        TDLen = FDLen - 1
        ReDim Dest(TDLen - 1)
        
        
        For n = 0 To TDLen - 1
            For k = 0 To FDLen - 1
                Dest(n) = Dest(n) + Src(k) * Cos(Pi / TDLen * n * k)
            Next k
        Next n
        Dest(0) = Dest(0) / 2
    End Sub
    While this works for all data sets, such that performing the DCT_B and then the IDCT_B on the output of the DCTB gives you back an exact copy of your original data (even if that data is a sinewave instead of a cosine wave), the ability to read the amplitude of a wave by looking at the value in the bin corresponding to its frequency works only for cosine waves. To be able to read off the amplitude of a sinewave, I have devised an equivalent DST, that I call the DST_B (Discrete Sine Transform, Ben's version) and its inverses, the IDST_B. Below is the code for the DST_B and IDST_B.

    Code:
    Public Sub DST_B(ByRef Src() As Double, ByRef Dest() As Double)
        Dim n As Long
        Dim k As Long
        Dim TDLen As Long
        Dim FDLen As Long
        
        TDLen = UBound(Src) + 1
        FDLen = TDLen + 1
        ReDim Dest(FDLen - 1)
        
        For k = 0 To FDLen - 1
            For n = 0 To TDLen - 1
                Dest(k) = Dest(k) + Src(n) * Sin(Pi / TDLen * n * k) / TDLen * 2
            Next n
        Next k
    End Sub
    
    
    
    Public Sub IDST_B(ByRef Src() As Double, ByRef Dest() As Double)
        Dim n As Long
        Dim k As Long
        Dim TDLen As Long
        Dim FDLen As Long
        
        FDLen = UBound(Src) + 1
        TDLen = FDLen - 1
        ReDim Dest(TDLen - 1)
        
        
        For n = 0 To TDLen - 1
            For k = 0 To FDLen - 1
                Dest(n) = Dest(n) + Src(k) * Sin(Pi / TDLen * n * k)
            Next k
        Next n
    End Sub
    The downside to using a DST_B is that the bin for DC is always 0, and so is the bin for the highest frequency. Interestingly, while the DST_B can produce frequency-domain data that allows reconstruction of a DC (despite the missing DC coefficient) or highest-frequency (despite the missing highest-frequency coefficient) signal via an IDST_B, the reconstructed time-domain signal always has its first sample set to 0.


    But what about measuing the amplitude of a sinusoidal wave form with an arbitrary phase shift (neither Sine nor Cosine)? For that, you simply take the bin corresponding to the frequency in question and take the DST_B copy of it and square it, and do the same thing for the DCT_B copy of it. Then you add these 2 squared values together and take the square root. This will give you the amplitude of any sinusoidal waveform with an arbitrary phase-shift. To make this process easier, I have created a Sub called DFT_B (which stands for Discrete Fourier Transform, Ben's version), which performs both the DCT_B and DST_B. I also created its inverse called the IDFT_B, which performs the IDCT_B and IDST_B and then averages those 2 results. Because both the DCT_B and DST_B both produce frequency-domain data that can reconstruct the entire signal (except for a DC signal when the transform is the DST_B), performing an IDFT_B isn't necessary, and I've just written it for the sake of completeness. Losing either the DCT_B or DST_B frequency-domain dataset only cuts the amplitude of the time-domain dataset by half when performing the IDFT_B, rather than changing the entire shape of the waveform.

    Here is my code for the DFT_B and IDFT_B.
    Code:
    Public Sub DFT_B(ByRef Src() As Double, ByRef DestDCT() As Double, ByRef DestDST() As Double)
        DCT_B Src(), DestDCT()
        DST_B Src(), DestDST()
    End Sub
    
    
    
    Public Sub IDFT_B(ByRef SrcDCT() As Double, ByRef SrcDST() As Double, ByRef Dest() As Double)
        Dim n As Long
        Dim DestPart2() As Double
        IDCT_B SrcDCT(), Dest()
        IDST_B SrcDST(), DestPart2()
        For n = 0 To UBound(Dest)
            Dest(n) = (Dest(n) + DestPart2(n)) / 2
        Next n
        Dest(0) = Dest(0) * 2
    End Sub
    Last edited by Ben321; Aug 10th, 2016 at 06:46 PM.

  2. #2

Posting Permissions

  • You may not post new threads
  • You may not post replies
  • You may not post attachments
  • You may not edit your posts
  •  



Click Here to Expand Forum to Full Width