Why does signal still contain a large number of spikes after FFT

I have perform FFT on my signal but due to the dc offset my amplitude peaks at 0 hz so i use x = x - mean(x) to remove the dc offset before fft but this time round after fft although I can see the difference but it still contains a lot of spikes. It is as if the signal has yet to go through FFT. Is there a reason for this ?

Answer 1

[dpb]

Well, if your original signal is full of a large number of frequency components, the spectrum will be, too. This looks like a few dominants with a large number of what may be sidebands and harmonics/subharmonics. You might see what a dB scale (log) looks like instead of just linear.

What is the signal? Knowing something of the system might let somebody here expound at depth; ya' never knows, but there's a "veritable plethora" of backgrounds/experience.

One issue that could be a problem depending -- was the data collected with analog anti-aliasing filters or some other way to ensure sampling was fast enough to avoid aliasing? If there were components in the original signal higher than or very near the Nyquist, they will fold back into the computed frequency range. If that is the case, once the data are sampled, there's no way to remove that source of contamination.

[crixus]

thanks for your reply. Actually my signal is measure of sound level from an arduino sensor given in decibels taken every second.

So i have several decibels readings around 1000-1500 in the range of maybe 70 + to 80+ decibels. I applied the fft to see what I can find but it turns out there are still a lot of spike after fft and removing the dc offset.

I tried do a log plot (loglog, log y, log x)

Is there any way to work around this ? or my signal is just like this ? Like using some form of filtering etc ?

[dpb]

"...is measure of sound level from an arduino sensor given in decibels taken every second."

Well, there's the issue as outlined before; anything above 0.5 Hz in the input is aliased back into the range. At 1Hz sampling rate, the Nyquist frequency is 0.5Hz. I don't know how you computed your frequency scale but it's incorrect; the Fmax is 0.5 Hz.

You don't say what sound you were recording, but again at the sampling frequency anything over 1 Hz is aliased so presuming it was speech or any thing in audible range, it's grossly undersampled and it's that effect that you're seeing in the spectrum.

Again as noted before, you have to sample at at least twice the rate of the highest frequency in the signal and it should be 3- or 4X. This can only be accomplished by either analog low-pass filtering prior to sampling or ensuring there's not any higher frequency components in the input signal. Once you've sampled it; there's no fixing it.

[dpb]

"... any way to work around this? or my signal is just like this? Like using some form of filtering, etc.?"

To reiterate, not having noticed the above plaint previously...

1. Fmax=1/2dt --> 1/2(1 sec) = 0.5 Hz (Nyquist frequency)
2. Therefore, anything in input signal over 0.5Hz is aliased
3. Once sampled, CANNOT undo the damage done
4. Ergo, this is "your signal" but it undoubtedly is useless to represent anything about whatever was the input to the sensor
5. Must sample at least 2X Fmax

[crixus]

Apologies, I'm a beginner in signal processing, having a hard time to understand the message you're conveying. I will look back at my code and try some things out before i ask anything

[dpb]

The code is NOT the issue; it's your data collection process.

Look up "Nyquist" and read the associated discussion you'll find there regarding sampling rates; there's where the problem lies. You simply cannot sample a 25Hz input at 1Hz and expect to get anything useful. See the example I posted in another comment.

[crixus]

Oh ok, because I have tried taken the readings at 0.1 s and even 0.01 second but the results are still as spiky which is well i thought the problem is with the code but I will go and try higher sampling rate to check the differences.

[crixus]

Hi dpb, thanks for your reply regarding the sampling rate. I look at my arduino code and realize no matter how fast i set the delay, the highest it can go is collecting data at every 0.3 - 0.4 second which it why the results of fft are the same. I'll go try to look up on increasing the sampling rate of my sensor. Thanks guys for your input

[dpb]

Also note that software-scheduled sampling is going to have "jitter" and that will grossly contaminate the results...the FFT assumes the sampling is precisely uniform dt and a little error in absolute time is likely a very large error in the actual sampling rate from one acquisition to the next.

To do this kind of thing you need an interrupt-driven system with very precise and tight handler code or, typically, one uses a dedicated DAQ to handle the data collection then does the processing of the sampled waveform.

I see in another response you've now recognized that simply setting a sample rate in software and doing an analysis of a set of data with that number has no bearing on what the real sample rate may have been...

[crixus]

Hi dpb, ok thanks for the tips, I greatly appreciate it. I will go try it out first and examine my output.

[dpb]

I'm not really at all familiar w/ Arduino (other than having heard of it) so looked at the documentation.

It appears it is only software-polled operation; altho what I saw says in a very tight loop it could approach 10 kHz doing nothing else.

Using it thru the Matlab serial port and with anything else going on will slow that down markedly plus the issue mentioned above of controlling "jitter" in the sampling rate that will, if it exists, wreak havoc on trying to do spectral analysis on the results.

I'm guessing if you could write a compiled routine that you dispatched to do nothing but collect a buffer of data, you might be able to make this work to a degree(*) to at least get a general feel for what the frequency response of the DUT is, but if there's a need for real accuracy in the results I think you'll need to move to some other DAQ designed for the purpose of collecting such data.

(*) If you were to go to the effort of writing sophisticated, interrupt-driven acq (assuming the Arduino supports the facility to be able to do that), up to whatever the actual hardware limits are, it's doable. When you start polling and running under an OS like Windows and the like that can and does swap processes at will, things are much tougher to time accurately enough for the purpose.

[Image Analyst]

Raspberry Pi and BeagleBoard are more sophisticated than Arduino (or so I've heard) so maybe you might need to switch to one of those.

[crixus]

Hi guys, thanks for the feedback ! truly appreciate it ! what i'm planning to do is to speed up the data saving portion and see what rate is can achieve, if there's a need i'll eliminate one of the board to lessen the data communicate between the microcontroller. otherwise I will just have to buy a real daq :( which i hope not (am a poor part time student)

[dpb]

Remember that for spectral analysis, besides the rate the consistency of that rate is important, too...sampling at 1 kHz a few fractions of a msec owing to OS task-swapping or garbage collection or whatever is a significant error in terms of the sample interval.

Answer 2

[Image Analyst]

What are you listening to? It looks like white noise. Is it? If so, I would not be surprised to see basically random energy in any of the wavelengths, which is what I see in your FFT. Like dpb said, we need to know what sound you expect your signal to be, and what you expect the noise to be, to best figure out what kind of filter you need to separate the two.

[crixus]

The sound that I am listening to is an underwater thruster running at approximately 1500rpm.

[Image Analyst]

Why are you surprised at the way it looks?

[dpb]

So the base frequency of the signal is 1500/60 --> 25 Hz. Ergo, you must sample at at least 50 Hz to catch that fundamental frequency and at sample rate if there were only a single frequency in the signal at that precise frequency you'd only be getting one value at each of two points; not much from which to interpolate a sine wave. Add any higher harmonics of which there are bound to be some and the actual frequency content will extend up to the mid-hundred Hz range or higher. Sampling at 1 second is nonsense.

Here's an example that numerically illustrates what happens when you do that...

>> dt=0.01/pi;  % a more-or-less randomized portion of time series
>> t0=0:dt:200; % from a non-integer dt...
>> s=sin(2*pi*1500/60*t0);  % and a signal at your dominant f
>> isamp=interp1(t0,1:length(t0),1:200,'nearest');  % now sample at 1 sec
>> isamp(end)=length(t0);  % clean up extrapolated point NaN from interp1
>> tsamp=t0(isamp);  % the times at 1 sec intervals
>> ssamp=s(isamp);   % associated data to the same times
>> Ssamp=fft(ssamp); % transform under-sampled wave
>> figure,plot(2*abs(Ssamp(1:length(Ssamp)/2+1))) % plot PSD
>>

Note that from a single pure, 25 Hz waveform, by undersampling we've gotten a PSD full of a bunch of peaks. The dominant is there, but so are a bunch of multiples that shouldn't be there at all; remember there was only a single pure tone in the actual signal.

All the rest are artifacts from the aliasing by picking only a point here and there from the true waveform with not enough resolution (short enough dt) to be able to draw a complete sine wave; if you plot ssamp from the above you'll see a roughly triangular waveform, not a sine.

As noted, once the signal has been sampled, there's absolutely no way to fix the problem. Digital filtering can do a lot to remove noise in the real signal or other artifacts but it cannot do anything at all to help in this case where the data are simply not at all representative of the actual waveform owing to (in this case) drastic undersampling.

[Image Analyst]

I'm not even sure what the "problem" is. He said "I have perform FFT on my signal" - well, for what purpose? OK, so he does FFT on that spiky signal and gets a spiky FFT. So, big deal? It is what it is. I'm not clear that the desired output is. Does he believe the FFT should not be what it is and want something smooth instead? I have no idea.

I also don't know exactly what "...is measure of sound level from an arduino sensor given in decibels taken every second." means? Does that mean you are taking one reading every second? Does that mean that every second you are taking thousands of readings and storing all of them in an array? I'm not sure. I'm hoping it's the latter because if you have an "underwater thruster running at approximately 1500rpm." I'm sure it's going to generate frequencies in the audible range so taking just one sample every second will not get you the spectrum of the sound as people would hear it. Like dpb says you'd only get super super low frequencies, way way lower than people could hear.

[dpb]

"Does that mean you are taking one reading every second? Does that mean that every second you are taking thousands of readings..."

I actually hadn't thought about of the latter, IA, but it is possible...if that were the case, then the likelihood of my supposition is markedly reduced, indeed.

And, whatever the "thruster" is, I suspect it is indeed likely to have a very rich output spectrum so OP's result may be what it actually is...not way to know certainly w/ the info given but I was presuming was just taking a sample once a second as initially written. I'd've thought the high-speed acquisition just downloaded once/sec. would be described differently; but who's to say???

[crixus]

Actually, this is my first time trying FFT and from the examples that I read, I thought i'm "supposed" to see only one or two spikes within a signal after FFT which is why i thought i'm wrong at this.

[dpb]

The code itself looks ok; it's the data-collection that is the issue as you've now recognized it seems. I'll simply note the "richness" of the FFT is simply dependent upon what the input content is; there are going to be as many peaks as there are discernible frequency components in the input signal; probably the examples in texts you've looked at were selected precisely to illustrate some principles from "clean" systems that had well-defined peaks of given frequencies.

In your particular case here, the system undoubtedly is very rich in energy content from the process and then by drastically under-sampling as your collection process did, that compounded the problem immensely by all of those various actual frequencies getting folded into the spectrum over a number of different times and so it ended up looking essentially as white noise as IA noted.

Answer 3

[crixus]

Hi both, thanks for your replies. The code that i used to perform the FFT is as follow

    x = xlsread('test.csv');
    x = x(~any(isnan(data1),2),:);
    x = x(1:1500,:);

    x = x-mean(x);

    fs=1000;                % sampling freq
    xmax = max(abs(x));     % find the maximum value
    x = x/xmax;             % scaling the signal so the largest value(xmax) is 1
    % time & discretisation parameters
    N = length(x);
    t = (0:N-1)/fs;      

    % spectral analysis
    win = hanning(N);       % sin window 
                            % default is symmetric, periodic start at 0 at
                            % point 1, periodic Hann window is constructed by extending 
                            % the desired window length by one sample, constructing a 
                            % symmetric window, and removing the last sample 
                            % symmetric starts at very close to zero ends at very close to zero
                            % hamming will being points slightly above 0, hanning starts at 0

    K = sum(win)/N;         % coherent amplification of the window

    X = abs(fft(x.*win))/N; % Without abs it will give complex number
                            % abs will give the complex magnitude
                            % FFT of the windowed signal

    NUP = ceil((N+1)/2);    % calculate the number of unique points
                            % ceil = round off to nearest positive infinity
    X = X(1:NUP);           % FFT is symmetric, throw away second half
    phase = unwrap(angle(fft(x.*win)));
    phase = phase(1:NUP);

    if rem(N, 2)            % odd nfft excludes Nyquist point
                            % rem = reminder after division
      X(2:end) = X(2:end)*2;
    else                    % even nfft includes Nyquist point
                            % The DC component and Nyquist component, 
                            % if it exists, are unique and should not be multiplied by 2
      X(2:end-1) = X(2:end-1)*2;
    end
    f = (0:NUP-1)*fs/N;     % frequency vector
    %X = 20*log10(X);        % spectrum magnitude

    % plotting of the spectrum
    figure()
    plot(f, X, 'r') % Create a plot with a logarithmic scale for the x-axis and a linear scale for the y-axis
    hold on
    xlim([0 max(f)]) % set the limit for x axis
    grid on
    set(gca, 'FontName', 'Times New Roman', 'FontSize', 14)
    title('Amplitude Spectrum of the Signal')
    xlabel('Frequency (Hertz)')
    ylabel('Magnitude')

    spectralcentroid = sum(X .* f')/sum(f);
    plot(spectralcentroid,'*','MarkerSize',14)
    hold off
    %%Plotting amplitude against phase
    figure()
    plot(f,phase/180*pi)
    set(gca, 'FontName', 'Times New Roman', 'FontSize', 14)
    title('Phase Spectrum of the Signal')
    xlabel('Frequency (Hertz)')
    ylabel('Phase (Radian)')

plotting of the spectrogram

    figure()
    len_seg=512; %window size 1024;
    spectrogram(x, len_seg, 3/4*len_seg, [], fs, 'yaxis')   % S = spectrogram(X,WINDOW,NOVERLAP,NFFT,Fs)
    h = colorbar;
    set(h, 'FontName', 'Times New Roman', 'FontSize', 14)
    ylabel(h, 'Magnitude (Decibels)');
    set(gca, 'FontName', 'Times New Roman', 'FontSize', 14)
    xlabel('Time (Seconds)')
    ylabel('Frequency (Hertz)')
    title('Spectrogram of the Signal')

[Image Analyst]

You forgot to attach "data" and the code to read it in, so no one can try your code.

[crixus]

oops, I have edited my post and attach the data

[dpb]

Doesn't make any difference; the problem is in the fundamental data collection and can't be cured except by sampling the input at (at least) 100 Hz, and probably preferably 2-3X times that.