Hi,

I've been trying out spectral analysis using FFT. I've modified one the sample programs given in MATLAB. Here's my code:

Fs = 100;                    % Sampling Frequency
T = 1/Fs;                    % Sample time
L = 400;                     % Length of signal
t = (0:L-1)*T;             % Time vector
x = (7+10*sin(2*pi*15*t));

NFFT = 2^nextpow2(L);     % ***PROBLEM***
fspec = Fs/2*linspace(0,1,NFFT/2+1);
X = fft(x,NFFT)/L;
Xf = 2*abs(X(1:NFFT/2+1));
Xf(1) = Xf(1)/2;                 % Restoring DC component to original value
stem(fspec,Xf);

The program works great! However if I try to increase NFFT, say:

NFFT = 4*2^nextpow2(L);

then the amplitudes in the immediate neighborhood of the DC component will be unusually high, distorting my beautiful spectrum: http://img841.imageshack.us/img841/2807/spectrum.png

These points might be due to the higher frequency resolution. My question is: How do I scale them properly?

I tried halving the first few components:

Xf(1:4) = Xf(1:4)/2;

but I don't think this a good solution. Advice?

On Aug 26, 4:54=A0pm, "Akhil Prem" <ramu.premku...@gmail.com> wrote:
> Hi,
>
> I've been trying out spectral analysis using FFT. I've modified one the s=
ample programs given in MATLAB. Here's my code:
>
> Fs =3D 100; =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0% Sampling Frequency
> T =3D 1/Fs; =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0% Sample time
> L =3D 400; =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 % Length of signal
> t =3D (0:L-1)*T; =A0 =A0 =A0 =A0 =A0 =A0 % Time vector
> x =3D (7+10*sin(2*pi*15*t));
>
> NFFT =3D 2^nextpow2(L); =A0 =A0 % ***PROBLEM***
> fspec =3D Fs/2*linspace(0,1,NFFT/2+1);
> X =3D fft(x,NFFT)/L;
> Xf =3D 2*abs(X(1:NFFT/2+1));
> Xf(1) =3D Xf(1)/2; =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 % Restoring DC compone=
nt to original value

The Nyquist component is zero. In general, however,
you also have to include

Xf(end) =3D Xf(end)/2;     % Restoring Nyquist component to original
value

> stem(fspec,Xf);
>
> The program works great! However if I try to increase NFFT, say:
>
> NFFT =3D 4*2^nextpow2(L);
>
> then the amplitudes in the immediate neighborhood of the DC component wil=
l be unusually high, distorting my beautiful spectrum:http://img841.imagesh=
ack.us/img841/2807/spectrum.png
>
> These points might be due to the higher frequency resolution.
> My question is: How do I scale them properly?
>
> I tried halving the first few components:
>
> Xf(1:4) =3D Xf(1:4)/2;
>
> but I don't think this a good solution. Advice?

the fft is linear

x1 =3D 7*ones(1,L);
x2 =3D 10*sin(2*pi*15*t));
for NFFT =3D [L, 2^nextpow2(L),  4*2^nextpow2(L) ]
    X1 =3D fft(x1,NFFT)/L;
    X2 =3D fft(x2,NFFT)/L;
    X  =3D X1+X2;
    etc % Compare the spectra of X1,X2 and X
end
Now, YOU decide what to do.

Hope this helps.

Greg

On Aug 26, 1:54=A0pm, "Akhil Prem" <ramu.premku...@gmail.com> wrote:
> Hi,
>
> I've been trying out spectral analysis using FFT. I've modified one the s=
ample programs given in MATLAB. Here's my code:
> ...

>
> The program works great! However if I try to increase NFFT, say:
>
> NFFT =3D 4*2^nextpow2(L);
>
> then the amplitudes in the immediate neighborhood of the DC component wil=
l be unusually high, distorting my beautiful spectrum:http://img841.imagesh=
ack.us/img841/2807/spectrum.png
>
> These points might be due to the higher frequency resolution. My question=
 is: How do I scale them properly?
> ...

These new points are scaled correctly. Your first L point fft just
didn't generate samples at those frequencies.

When you generate a finite set of (L) samples of a sum of tones, the
frequency response is the frequency response of the tones convolved
with the frequency response of the (in this case rectangular) window
which is a sinc function. When the window is L points long, you
calculate an L point fft, and the tones are bin centered in that fft,
the sync functions have zeros at all output samples except the
positions of the tones. Calculate larger and larger ffts and you
sample more and more non-zero points on the frequency response of the
windowed sum of tones. If you try L points and 2*L points you will
generate, first all the possible zero samples and non-zero samples at
tones and second all those samples plus all the non-zero samples
halfway in between the first samples.

> Advice?

Your code is correct, learn to understand what it means.

Dale B. Dalrymple

dbd <dbd@ieee.org> wrote in message <d595b3c8-ef90-4d3b-b85d-0a01113b1728@p12g2000prn.googlegroups.com>...
> On Aug 26, 1:54 pm, "Akhil Prem" <ramu.premku...@gmail.com> wrote:
 
> These new points are scaled correctly. Your first L point fft just
> didn't generate samples at those frequencies.
> 
> When you generate a finite set of (L) samples of a sum of tones, the
> frequency response is the frequency response of the tones convolved
> with the frequency response of the (in this case rectangular) window
> which is a sinc function. When the window is L points long, you
> calculate an L point fft, and the tones are bin centered in that fft,
> the sync functions have zeros at all output samples except the
> positions of the tones. Calculate larger and larger ffts and you
> sample more and more non-zero points on the frequency response of the
> windowed sum of tones. If you try L points and 2*L points you will
> generate, first all the possible zero samples and non-zero samples at
> tones and second all those samples plus all the non-zero samples
> halfway in between the first samples.
> 
> > Advice?
> 
> Your code is correct, learn to understand what it means.
> 
> Dale B. Dalrymple

Thank you Greg and Dale.

First of all, I'm aware of why this is happening - about sinc functions and windows, as you have rightly explained. What I'm looking for is a solid technique to handle these unruly samples despite any value of NFFT I may choose. In other words, is there a relation between frequency resolution and the number of samples to be scaled (halved?) near the DC value.

Greg's solution is interesting, but shouldn't I have knowledge of the DC value beforehand? Maybe mean() would do...

I've fixed it.

I took the fft of the mean of the given signal as Gereg suggested and subtracted it from my fault spectrum. He're the code:

Fs = 100;
T = 1/Fs;                     % Sample time
L = 400;                      % Length of signal
t = (0:L-1)*T;              % Time vector
x = (7+10*sin(2*pi*15*t));
xdc = mean(x)*ones(1,L);

NFFT = 2*2^nextpow2(L);
fspec = Fs/2*linspace(0,1,NFFT/2+1);
X = fft(x,NFFT)/L;
Xdc = fft(xdc,NFFT)/L;
Xf = abs(2*(X(1:NFFT/2+1))-(Xdc(1:NFFT/2+1)));  % FIXED

Thanks for the help.

On Aug 26, 11:44=A0pm, "Akhil Prem" <ramu.premku...@gmail.com> wrote:
> ...
>
> Thank you Greg and Dale.
>
> First of all, I'm aware of why this is happening - about sinc functions a=
nd windows, as you have rightly explained.

If this sentence were true, you wouldn't state the next two sentences.

> What I'm looking for is a solid technique to handle these unruly samples =
despite any value of NFFT I may choose.

Try plotting only the values that are peaks of the periodogram or use
"interpolated fft" or frequency reassignment" algorithms to get more
accurate peak locations than peak picking gives for non-bin-centered
components.

> In other words, is there a relation between frequency resolution and the =
number of samples to be scaled (halved?) near the DC value.
>

Only the DC and Nyquist samples (1 each) may deserve special scaling.

> Greg's solution is interesting, but shouldn't I have knowledge of the DC =
value beforehand? Maybe mean() would do...

Since your spectral components are all bin centered, the DC component
and the mean are the same in either the time or frequency domains.

Dale B. Dalrymple

> Only the DC and Nyquist samples (1 each) may deserve special scaling.

This is where the problem arises. Its NOT just the DC or Nyquist component that needs special scaling if one uses normal fft (without ""interpolated fft" or frequency reassignment" algorithms"). A few components surrounding the DC value need also be scaled. My question was, how do I figure out which all components I should scale.

Anyway, I've found a solution I'm satisfied with. Thanks again.

On Aug 27, 9:59=A0am, "Akhil Prem" <ramu.premku...@gmail.com> wrote:
> > Only the DC and Nyquist samples (1 each) may deserve special scaling.
>
> This is where the problem arises. Its NOT just the DC or Nyquist componen=
t that needs special scaling if one uses normal fft (without ""interpolated=
 fft" or frequency reassignment" algorithms"). A few components surrounding=
 the DC value need also be scaled. My question was, how do I figure out whi=
ch all components I should scale.

I think you have found the best solution for yourself. When you have
added to your artificial signal a component that annoys or confuses
you, the best way to remove the annoyance and confusion is to subtract
that signal component out. If you move from artificial data to real
world sampled data you may need different solutions. Good luck!

>
> Anyway, I've found a solution I'm satisfied with. Thanks again.

The source of your bliss is obvious.

Dale B. Dalrymple

> If you move from artificial data to real
> world sampled data you may need different solutions. Good luck!

This program actually forms a part of a home brewed data acquisition system I'm constructing and it is giving good results with some real world signals I've sampled:

http://img266.imageshack.us/img266/3600/spectrumof555waveforms.png

Can I ask for more?

On Aug 27, 12:59=A0pm, "Akhil Prem" <ramu.premku...@gmail.com> wrote:
> > Only the DC and Nyquist samples (1 each) may deserve special scaling.
>
> This is where the problem arises. Its NOT just the DC or Nyquist componen=
t that needs special scaling if one uses normal fft (without ""interpolated=
 fft" or frequency reassignment" algorithms"). A few components surrounding=
 the DC value need also be scaled. My question was, how do I figure out whi=
ch all components I should scale.
>
> Anyway, I've found a solution I'm satisfied with. Thanks again.

I'm not sure what your application is. However, the fft,
with appropriate windowing, reveals the frequency content
of your signal. If you fiddle with the spectrum you might get
something that looks pleasing to the eye but the corresponding
frequency content is just a figment of your imagination.

I you don't want sidelobes of a large DC component to
contaminate the rest of the spectrum, subtract the mean
from the time signal instead of fiddling with the spectrum
a posteriori.

Hope this helps.

Greg

> I'm not sure what your application is. However, the fft,
> with appropriate windowing, reveals the frequency content
> of your signal. If you fiddle with the spectrum you might get
> something that looks pleasing to the eye but the corresponding
> frequency content is just a figment of your imagination.

Hmm... you and dale may be right. The problem is that I don't have a standard spectroscope or the  sort to compare my results with. I'll have to investigate more.

> I you don't want sidelobes of a large DC component to
> contaminate the rest of the spectrum, subtract the mean
> from the time signal instead of fiddling with the spectrum
> a posteriori.

But by linearity property, both approaches should yield the same results, shouldn't it? And I do wish to see the DC value in my spectrum also.

Here's another test I conducted; this time I sampled waveforms from an accelerometer. The waveforms are shaking it:

http://img839.imageshack.us/img839/9985/accelerometerxandychann.png

Ultimately, my application boils down to this: I want to select an optimal sampling frequency for sampling an accelerometer signal and I'm trying to analyse its spectrum to that end.

On Aug 28, 10:11=A0am, "Akhil Prem" <ramu.premku...@gmail.com> wrote:
> > I'm not sure what your application is. However, the fft,
> > with appropriate windowing, reveals the frequency content
> > of your signal. If you fiddle with the spectrum you might get
> > something that looks pleasing to the eye but the corresponding
> > frequency content is just a figment of your imagination.
>
> Hmm... you and dale may be right. The problem is that I don't have a stan=
dard spectroscope or the =A0sort to compare my results with. I'll have to i=
nvestigate more.
>
> > I you don't want sidelobes of a large DC component to
> > contaminate the rest of the spectrum, subtract the mean
> > from the time signal instead of fiddling with the spectrum
> > a posteriori.
>
> But by linearity property, both approaches should yield the same results,=
 shouldn't it? And I do wish to see the DC value in my spectrum also.
>
> Here's another test I conducted; this time I sampled waveforms from an ac=
celerometer. The waveforms are shaking it:
>
> http://img839.imageshack.us/img839/9985/accelerometerxandychann.png
>
> Ultimately, my application boils down to this: I want to select an optima=
l sampling frequency for sampling an accelerometer signal and I'm trying to=
 analyse its spectrum to that end.

Concentrate on

X2 =3D fft(x-mean(x),NFFT)/length(x);

Hope this helps.

Greg