Construct a sum of N trigonometric functions to reach a given global period

Question

I have programmed a model based on a research article. In order to validate it I need to recreate the same dataset as in the article.

My goal is to construct a signal by superposition of a large number, $N$, of sinusoidal functions with periods ranging between $T_{min}$ and $T_{max}$, with a constant amplitude and random phase (bandpass filter). The frequencies of these functions would be "chosen" in order to reach for the summed signal a given period $T_\Sigma$, way larger than $T_{max}$.

I have no problem to deduce from a sum of sinusoidal functions what would be the period of the signal : $T_\Sigma = lcm\left[T_1;T_i;T_N\right]$, where lcm is the least common multiple of every individual period.

But doing the opposite: building a list of numbers to achieve a given lcm seems harder to me...

To sum everything up, my signal would look like this:

$$\sum_{i=1}^N \sin\left(2\pi\frac{1}{T_i}t\right)$$ Where $T_{min}<T_i<T_{max}$

Thus the period of the summed signal would be $T_\Sigma = lcm\left[T_1;T_i;T_N\right]$.

Fixed value being $T_{min}$, $T_{max}$, $T_\Sigma$ and $N$, I need to find a relationship allowing me to pick every $T_i$.

Thank You !

Sidenote: yes, I know I could generate a bandpass signal based on filtered white noise, but then it would not have any "global" period, as it was performed in the article results I'm trying to reproduce.

Answer 1

Let the prime factorization of $T_\Sigma$ be $ p_1^{\alpha_1} p_2^{\alpha_2} \cdots p_q^{\alpha_q}. $

For each $p_i^{\alpha_i}$, you need at least one $T_j$ that is a positive multiple of $p_i^{\alpha_i}$. The positive multiple could be $p_i^{\alpha_i}$ itself, but if you never have a $T_j$ with at least $\alpha_i$ factors of $p_i$ then the LCM will not contain $\alpha_i$ factors of $p_i$ and will not be $T_\Sigma$.

A single $T_i$ could be a multiple of some $p_i^{\alpha_i}p_j^{\alpha_j}$, that is, it could cover two of the primes at once, so you do not necessarily need a different $T_i$ for every prime in the prime factorization. This is how you can deal with the case where $q > N,$ that is, if you have more distinct primes than distinct periods.

Note that if $T_\Sigma$ is a power of a prime, for example, if $T_\Sigma = 1024,$ the only way to achieve that period is for $T_\Sigma$ itself to be one of the periods in your list.

To satisfy the constraint that $T_\min \leq T_i \leq T_\max,$ you may have to multiply some of the $p_i^{\alpha_i}$ by other factors from the prime factorization in order to get a period between $T_\min$ and $T_\max.$ If $q > N$ you may have to try different combinations of $p_i^{\alpha_i}p_j^{\alpha_j}$ with other factors in order to find one that "fits." Note that if $p_i^{\alpha_i} > T_\max$ for any $i$ you are out of luck. But even if that does not happen, a particular problem still may be unsolvable for other reasons.

One approach might be to make a list of all factors of $T_\Sigma$ that are between $T_\min$ and $T_\max.$ There are only $\prod_i (\alpha_i + 1)$ factors, which you can generate systematically and discard any that are out of the desired range. You can also (for now) set aside any factor that is not a multiple of $p_i^{\alpha_i}$ for at least one $p_i.$ Check the list of factors to make sure each $p_i^{\alpha_i}$ divides at least one member of the list. Now choose a factor from the list that is a multiple of $p_1^{\alpha_1},$ a factor that is a multiple of $p_2^{\alpha_2},$ and so forth, up to $N$ factors, until all the distinct prime factors are accounted for. (If $q > N$ you will need to observe when a chosen factor accounts for multiple distinct prime factors or you will not succeed in accounting for all of them.)

It may help if you make a memo of what factors of $p_i^{\alpha_i}$ each number in the list has. But even with aids like that, if $q > N$ you might have to backtrack (discard the choices made so far and try new ones) in order to find factors that account for every distinct prime factor without using more than $N$ factors of $T_\Sigma.$

If you are able to construct a list of $N$ or fewer factors that account for every distinct prime factor of $T_\Sigma,$ you can fill out the list with any other factors of $T_\Sigma$ between $T_\min$ and $T_\max,$ whether or not each new factor is a multiple of some $p_i^{\alpha_i}$.

Answer 2

The answer of @David K seems mathematically correct but didn't really suited my needs.

Instead I decided to built from scratch a DFT. I first generated a bandpass filter with a zero phase for every frequency sampled. Then I randomized the phase, and finally, I rebuilt my signal with a iDFT.

To be honest, I have no idea how to prove that this is mathematically correct, except maybe for one part: I'm using the periodicity property of the DFT to ensure the period of my final signal $T_\Sigma$. But I think it somehow respect the properties enounced by @David K

Since I'm far from being an expert and to be sure I don't use wrong general notation or definition, here is the MATLAB code I created to generate my signal.

%%% Signal parameters

    % Lowest frequency
    fL = 1;
    % Highest frequency
    fH = 2;
    % Sampling frequency
    Fs = 100;   
    % Sampling period
    T = 1/Fs;         
    % Total time of signal
    TSigma = 50;
    % Length of signal
    L = TSigma/T; 
    % Time vector
    time = (0:L-1)*T;  
    % Frequencies vector
    f = Fs*(0:(L/2))/L;

%%% Constructing the 'real' FFT (zero phase on the whole frequencies domain)

    % Initialization

    P1 = zeros(1,L/2+1);
    P2 = zeros(1,L);

    % Building the single-sided spectrum P1 with a real signal

    for i = 1:length(f)
        if(f(i)>= fL & f(i)<= fH)
            P1(i) = 1;
        end
    end

    P1(2:end-1) = P1(2:end-1)/2;

    % Building the two-sided spectrum P2

    P2(1:L/2+1) = P1;
    P2(L/2:L) = flip(P1,2);

    % Building the fft of the signal
    signalFFT = P2*L;

%%% Randomizing the phase

    % Creating a random vector with values between [-pi;pi]
    rndTheta= -pi + (2*pi).*rand(1,L/2-1); 

    % Adding a random phase noise to the signal one one side of the
    % spectrum
    signalFFT(:,2:L/2)=signalFFT(:,2:L/2).*exp(1i*rndTheta);

    % Adding the conjuguate of these values on the other side of the
    % spectrum
    signalFFT(:,L/2+2:L)=signalFFT(:,L/2+2:L).*exp(-1i*flip(rndTheta,2));

%%% Rebuilding the signal with computing the iDFT of the previouly built DFT 

    signal = ifft(signalFFT,'symmetric');

So I don't really know what I'm doing, but I'm pretty sure this is mathematically right.. If anyone want to help the community and prove/correct it, I would be very thankful :)

Answer 3

If I understand correctly what you asked, you can choose the $f_{i} = \frac{1}{T_{i}}$ from the list:

$$f_{i} = f_{min} + if_{\Sigma}, $$

for $i = 0,1,...,\frac{f_{max}-f_{min}}{f_{\Sigma}}$. Choose $i = 0, 1, 2,..., N-2$ and the final frequency choose $f_{max}$.