Unexpected peaks of FFT of positive-only sparse signal [OCTAVE or MATLAB]

Question

Suppose a signal corresponds to day values over a year (365 days). It is composed of all zeros except for a few sparse values which corresponds to isolated peaks all separated by the same interval (30 days). I obtain the frequency spectrum with the Fast Fourier Transform function.

1. how to get rid of the high 0Hz peak? EDIT: this is related to the non-zero-mean nature of the signal. See this post for more details.

2. The first peak is then at 12Hz, which is somehow expected. However, peaks are also present at 24Hz, 36Hz, 48Hz ... . Is it an aliasing problem? How to get rid of it?

Below is my code. It was tested in Octave but it should also work in Matlab

close all
clear all

T = 1/365; % period
samp_freq = 1/T;  % sample frequency
t=0:T:2; % overall time span is two years

% build signal
x= zeros(length(t),1);    
for i=1:length(t)
    if mod(i,30) == 0
       x(i) = 100; 
    else
        x(i) = 0;
    end
end

figure(1)
plot(t,x)
grid
xlabel("Time [years]")
ylabel("Signal amplitude")


y=fft(x);
N = length(x);

for i=1:N
    f(i) = (i-1)*samp_freq/N;
end

figure(2)
plot(f,abs(y))
xlabel("Frequency")
ylabel("Signal amplitude")

figure(3)
plot(f(1:80),abs(y(1:80)))
xlabel("Frequency")
ylabel("Signal amplitude")

Answer 1

The peak at 0 Hz represents the "DC component" of your signal, which is the same as the arithmetic mean. You could subtract the mean value from your data set to get rid of this, or just ignore it.

The other peaks at multiples of 12 Hz are called "harmonics" - they exist because your periodic signal is not sinusoidal (it's actually a pulse train). The first peak, at 12 Hz, is called the "fundamental". You could either pre-filter your signal to remove the harmonics, or just ignore everything above, say, 20 Hz).

You might also want to consider applying a window function to your data, in order to reduce spectral leakage and thereby make your peaks sharper and more accurate.

Answer 2

Your results are completely correct.

Your input signal is basically a dirac comb, with a period T of 30/(365*2): 100 III_{30/(365*2)}.

The exact fourier transform of this signal is again a dirac comb: 100 * (365*2)/30 III_{(365*2)/30} = 24333 III_{(365*2)/30}. This is almost the result that you get.

If you want to search the period of your input signal, search for the first peak after 0Hz.