5
$\begingroup$

I have been posting all over SE looking for a good method of finding the maximum of a trigonometric polynomial. A trigonometric polynomial is the sum of a set of sinusoids,

$$a_0 + \sum_{k = 1}^{N} a_k \cos(k \theta + \phi_k)$$

and can also be written in the form

$$ \sum_{k=-N}^{N} c_k \, e^{i k \theta} $$ where $c_k$ = $\bar{c}_{-k}$

Substituting $z = e^{i \theta}$ turns this into the typical form of a complex valued polynomial.

The method I am currently using is to take the find the roots of the derivative of the polynomial, then evaluate the original polynomial at these locations to find the maximum. While quite fast with MATLABs roots command, it is not quite fast enough.

I need a method where I can trade accuracy for speed. Does anyone know of any other fast methods? This is a similar problem to estimating the crest factor of a signal. In my problem $N$ ranges from 3 to 15.

Improve this question
$\endgroup$
5
  • 2
    $\begingroup$ Your question is a special case of a problem that has been studied quite extensively in the last 30 years or so. Search for crest factor, peak-to-average ratio (PAR), and peak to average power ratio (PAPR) for methods that might be used. $\endgroup$
    – Dilip Sarwate
    Commented Mar 20, 2013 at 21:32
  • $\begingroup$ What is $n$ in the equation you gave above? I'm guessing that the sinusoids aren't all at the same frequency. Should $n$ be $k$ instead? $\endgroup$
    – Jason R
    Commented Mar 20, 2013 at 22:31
  • $\begingroup$ @JasonR yep I have corrected. $\endgroup$
    – geometrikal
    Commented Mar 22, 2013 at 17:05
  • $\begingroup$ @DilipSarwate I have some papers on crest factor but haven't found any with good method, do you know one? $\endgroup$
    – geometrikal
    Commented Mar 22, 2013 at 17:05
  • $\begingroup$ @geometrikal As my comment above hinted, there are no good methods (that will work well at finding the crest factor of all signals). All known methods are flawed in some way or other. In each case, one can find examples where the chosen method works poorly, and also examples where the method works reasonably, quite well, or superbly. $\endgroup$
    – Dilip Sarwate
    Commented Mar 22, 2013 at 21:48

2 Answers 2

Reset to default
1
$\begingroup$

I would proceed as follows:

1) Implement trigonometric polynomial at equidistant 'theta' nodes by means of an 'fft'. Sampling density on 'theta' axis can be increased by interleaving zeros between subsequent elements of the input vector e.q c[-N] 0 0 0 c[-N+1] 0 0 0 ... c[N-1] 0 0 0 c[N]

2) Call 'max' function of Matlab to find maximal value and its index.

This should compute reasonably fast and you can trade the speed for accuracy by increasing the number of zeros.

UPDATE:

I have to correct myself - sampling density can be increased by zero-padding vector $c$ rather than interleaving zeros.

Following Matlab script demonstrates the idea by computing approximated values of a trigonometric polynomial by means of FFT.

% Create 'c' coefficients of a trigonometric polynomial approximation of a
% shifted sinc function
N = 6;
shift = -pi/N;
cc =  fftshift(ifft(conj(sinc(linspace(-2*pi+shift,2*pi+shift,2*N+1)))));

% Evaluate trigonometric polynomial computing it directly according to definition
% Maximal value of this trigonometric polynomial is ~0.923
theta = [0:-0.01:-2*pi];
[NN,THETA]=meshgrid([-N:N], theta);
y = (exp(i.*NN.*THETA)*cc.').';

% Plot the densely sampled trigonometric polynomial
figure;
plot(theta,real(y));

% Sample trigonometric polynomial at equidistant nodes using FFT
% theta = [-2*pi*[0:N-1]./N]
fft_approx{1} = conj(fft(ifftshift(cc)));

% Overlay the plot
hold on;
plot(-2*pi*[0:2*N]./(2*N+1),fft_approx{1},'.:r');

% Find a maximum value 
max_val{1} = max(real(fft_approx{1}));

% Increase the sampling density 5 times using zero-padding
os_factor = 5;
zpad = zeros(1,numel(cc)*floor(os_factor/2));
cc_padded = [ zpad cc zpad];

% Sample trigonometric polynomial at a denser grid using FFT
fft_approx{2} = conj(fft(ifftshift(cc_padded)));

% Overlay the plot
hold on;
plot(-2*pi*[0:numel(fft_approx{2})-1]./numel(fft_approx{2}),real(fft_approx{2}),'.:g');

% Find a maximum value 
max_val{2} = max(real(fft_approx{2}));

% Annotate the graph
legend('Original trigonometric polynomial',...
       sprintf('Sampled at %d points using fft',numel(fft_approx{1})),...
       sprintf('Sampled at %d points using fft',numel(fft_approx{2})));

% Display approximated maximal values   
fprintf('-----------------------------------------------------\n');   
fprintf('%d-points fft approximation of the maximal value = %g\n',numel(fft_approx{1}),max_val{1});
fprintf('%d-points fft approximation of the maximal value = %g\n',numel(fft_approx{2}),max_val{2});

fprintf('-----------------------------------------------------\n');
fprintf('True maximal value to 3 decimal places is 0.923\n');
Improve this answer
$\endgroup$
7
  • 2
    $\begingroup$ I don't understand how your method works, and request that you post an illustration to show what it does. To my mind, the FFT will give the amplitudes of the individual sinusoids (the $a_i$ and the $\phi_i$ are known already), whereas what the OP is looking for is the maximum value of the sum of those sinusoids. $\endgroup$
    – Dilip Sarwate
    Commented Mar 21, 2013 at 18:28
  • $\begingroup$ @Dilip DFT is by definition a trigonometric polynomial sampled at equidistant nodes. What I propose is to sample values of $$ \sum_{k=-N}^{N} c_k \, e^{i k \theta} $$ at nodes $$ -2\pi*\frac{[0:2N]}{2N+1} $$ using FFT (which is fast) and then find the maximum. $\endgroup$
    – Ilya Sedelnikov
    Commented Mar 22, 2013 at 12:59
  • $\begingroup$ @Dilip I also updated my earlier reply with a matlab script that implements the idea $\endgroup$
    – Ilya Sedelnikov
    Commented Mar 22, 2013 at 13:16
  • $\begingroup$ @IlyaSedelnikov I have compared your using FFT and simply evaluating the poly at the same points. There doesn't seem to be much difference in speed. For example, evaluating poly at 13 points (for loop summing c_k exp(i*k*(theta+t)) where t = nodes above) takes 0.000090 seconds, using FFT (conj(fft(ifftshift(cc))) from your code) takes 0.000089 seconds. Is DFT usually faster that plain matlab for loop / sum in your experience? $\endgroup$
    – geometrikal
    Commented Mar 22, 2013 at 17:38
  • $\begingroup$ @geometrikal: DFT is faster than direct evaluation as the number of sampling points grows. The number of points where the polynomial should be sampled depends on your target precision. It is probably worth comparing runtime of a 'roots' command on your typical inputs (N between 5 and 13) to the runtime of 'fft' on a number of points that gives you desired precision. For example, on my computer 'fft' applied to 130-element vector is twice faster than 'roots' applied to 13-element vector, but I do not know whether resulting precision will be acceptable for your purposes. $\endgroup$
    – Ilya Sedelnikov
    Commented Mar 22, 2013 at 18:05
1
$\begingroup$

This function is periodic with 2*pi and it it's highest frequency component is N times the fundamental frequency. As such it can be sampled at a rate of 1/2*N without losing information. However, there is no guarantee that sampling will hit the maximum or that maximum sampled value is close to the analog maximum. The likelihood of this increases drastically by oversampling, e.g sample 4*N or 8*N values.

So one thing to try would be to basically calculate the function values on a equidistant grid of, say 8*N, values, locate the maximum and then do a binary search around the maximum until the desired accuracy is reached.

Again there is no guarantee that this will ALWAYS work, however you can try it with a few actual data sets and see how much over sampling is required and whether this is actually faster than the root search.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.