View Raw SPL
/*****************************************************************************
* *
* AMPSPEC.SPL Copyright (C) 2010 DSP Development Corporation *
* All Rights Reserved *
* *
* Author: Randy Race *
* *
* Synopsis: N point normalized complex amplitude spectrum *
* *
* Revisions: 18 Jul 2010 RRR Creation *
* *
*****************************************************************************/
#if @HELP_AMPSPEC
AMPSPEC
Purpose: Calculates an N point complex amplitude spectrum.
Syntax: AMPSPEC(s, N, type)
s - The input series or array.
N - An optional integer, the amplitude spectrum length.
Defaults to the length of the input series.
type - Optional. A string, the output type.
"single" : single sided display (default)
"double" : double sided display
"shift" : double sided display shifted about 0 Hz.
Returns: A complex series or array, the N point normalized complex
spectrum of the input.
Example:
W1: gcos(1000, 1/1000, 100)
W2: ampspec(w1)
W3: ampspec(w1, "double")
W4: ampspec(w1, "shift")
W2 contains 500 complex values with a peak at 100 Hz. W3 contains
1000 values with peaks at 100 Hz and 900 Hz. W4 contains 1000
values with peaks at -100 Hz and +100 Hz. In all cases, the
magnitude of the peak values is 0.5, 1/2 the amplitude of the
input cosine.
calculated, the computation is performed more quickly.
Example:
W1: gsqr(1000, 1/1000, 4)
W2: ampspec(w1, "shift")
W3: {mean(mag(w1)^2)}
W4: {sum(mag(w2)^2)}
W3 == W4 == 0.5 verifying a form of Paresval's Theorem.
Remarks:
AMPSPEC computes N equally spaced samples of the normalized
complex amplitude spectrum by using the FFT. The raw FFT
values are normalized by the length of the input series
such that:
ampspec(s) = fft(s) / length(s)
For a sampling rate Fs, the default single sided amplitude
spectrum displays N/2 frequency values from 0 to Fs/2. The
double sided amplitude spectrum, "double", displays N values
from 0 to Fs and the shifted spectrum, "shift", displays
N values from -Fs/2 to Fs/2.
Unlike SPECTRUM, AMPSPEC does not scale the values between
0 and Fs/2 by 2. Thus, the single sided amplitude
spectrum of a 1 volt sinusoid of frequency F shows a peak
of 0.5 at frequency F and the double sided shifted amplitude
spectrum shows peaks of 0.5 located at -F and +F.
The "double" sided or "shift" amplitude spectrum obeys a
form of Parseval's Theorem such that:
mean(mag(s)^2) == sum(mag(ampspec(s)^2)
See MAGSPEC to display the magnitude spectrum.
See PHASESPEC to display the phase spectrum.
See SPECTRUM to compute a normalized frequency spectrum such
that a 1 volt sinusoid at frequency F displays a peak of
1 at frequency F.
See Also:
fft
magspec
phasespec
spectrum
#endif
/* returns complex amplitude spectrum - FFT(s) / length(s) */
ITERATE ampspec(s, N, type, func)
{
local a, L;
if (isunspecified(func))
{
func = __FUNC__;
}
(s, N, type) = ampspec_parse_args(s, N, type, func);
L = length(s);
if (N < 0) N = L;
/* normalized FFT */
a = fft(s, N) / L;
if (type == "single")
{
/* single sided spectrum */
a = extract(a, 1, int(N/2) + 1);
}
else if (type == "shift")
{
/* double sided - shift about DC */
a = fftshift(a);
}
return(a);
}
/* parse args */
ampspec_parse_args(s, N, type, func)
{
local temp;
if (argc < 3)
{
if (argc < 2)
{
if (argc < 1) error("ampspec - input series required");
N = -1;
}
type = "single";
}
if (isstring(N))
{
temp = N;
N = isstring(type) ? -1 : type;
type = temp;
}
type = ampspec_parse_string(type, func);
return(s, N, type);
}
/* parse string options */
ampspec_parse_string(str, func)
{
local frange = "single";
if (isstring(str))
{
if (strlen(str) > 0)
{
str = tolower(str);
switch (str)
{
case "onesided":
case "half":
case "single":
frange = "single";
break;
case "twosided":
case "whole":
case "double":
case "full":
frange = "double";
break;
case "center":
case "centered":
case "centerdc":
case "shift":
frange = "shift";
break;
default:
error(sprintf("%s - unknown flag '%s'", func, str));
}
}
}
return(frange);
}