MATLAB: Fourier Transformation of dirac comb

dirac combfftfourier transformationimpuls functionsignal processing

Hello guys,

in my code I generated a dirac comb and its FFT with:

%Time Signal
CarrierFrequenz=100; %frequency of the impulse in Hz
fs=CarrierFrequenz*10; % sample frequency (10 times higher than Carrier Frequency)
L=fs; %length signal
t=0:1/fs:1; % time vector
y=zeros(size(t));
y(1:fs/CarrierFrequenz:end)=1; %dirac comb
subplot(3,1,1)
plot(t,y); %time signal
xlabel('time(s)')
ylabel('P1')
%Frequency signal
Y = fft(y);
P2 = abs(Y/L);
subplot(2,1,2)
f = (fs*(-L/2:((L/2)-1))/L); %frequency domain
plot(f,P2(1:end-1)) 
title(' Amplitude Spectrum of S(t) FFT')
xlabel('f(Hz)')
ylabel('|P1(f)|')

and it generates the following plot:

Signal processing says:

So I wonder why the amplitudes in the frequency spectrum have this stair form and their magnitudes are way less than in the time domain.

Would be more than happy if s.o. could help me with that.

Best regards

Johannes

Best Answer

There is one too many points in the signal y. If you

plot([y y])

and zoom in you'll see the pulse at sample 1001 immediately followed by a pulse at sample 1002. So the underlying infinte length sequence isn't the infinite length sequnece in your question. But if you make y one sample shorter, then the underlying infinite sequence is the sequence in your question and it yields the expected result, at least with respect to eliminating the spreading (spectral leakaage) and non-constant peaks in the frequency domain. I didn't look carefully at the scaling.

CarrierFrequenz=100; %frequency of the impulse in Hz
fs=CarrierFrequenz*10; % sample frequency (10 times higher than Carrier Frequency)
L=fs; %length signal
t=0:1/fs:1; % time vector
t = t(1:end-1);  % this line changed

y=zeros(size(t));
y(1:fs/CarrierFrequenz:end)=1; %dirac comb
subplot(3,1,1)
plot(t,y); %time signal
xlabel('time(s)')
ylabel('P1')
%Frequency signal
Y = fft(y);
P2 = abs(Y/L);
subplot(2,1,2)
f = (fs*(-L/2:((L/2)-1))/L); %frequency domain
plot(f,P2(1:end)) % this line changed
title(' Amplitude Spectrum of S(t) FFT')
xlabel('f(Hz)')
ylabel('|P1(f)|')

Related Solutions

MATLAB: Time domain and frequency domain signal generation and plot

hello

here you are

hope it helps

clc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%





% load signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% data
% data = importdata('SFC5S_nov25_ST_1_1.txt');
data = importdata('Test_1.csv');
dt = 1/8000; % 

signal = data(:,1);
samples = length(signal);
Fs = 1/dt; % sampling frequency (Hz)
%% decimate (if needed)
% NB : decim = 1 will do nothing (output = input)
decim = 1;
if decim>1
    signal = decimate(signal,decim);
    Fs = Fs/decim;
    samples = samples/decim;
    
end
time = (0:samples-1)*1/Fs;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = 64;    % 
OVERLAP = 0.95;
% spectrogram dB scale
spectrogram_dB_scale = 80;  % dB range scale (means , the lowest displayed level is XX dB below the max level)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% options 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% if you are dealing with acoustics, you may wish to have A weighted
% spectrums 
% option_w = 0 : linear spectrum (no weighting dB (L) )
% option_w = 1 : A weighted spectrum (dB (A) )
option_w = 0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



% display 1 : time domain plot
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1),plot(time,signal,'b');grid
title(['Time plot  / Fs = ' num2str(Fs) ' Hz ']);
xlabel('Time (s)');ylabel('Amplitude');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 2 : averaged FFT spectrum
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[freq, sensor_spectrum] = myfft_peak(signal,Fs,NFFT,OVERLAP);
% convert to dB scale (ref = 1)
sensor_spectrum_dB = 20*log10(sensor_spectrum);
% apply A weigthing if needed

if option_w == 1
    pondA_dB = pondA_function(freq);
    sensor_spectrum_dB = sensor_spectrum_dB+pondA_dB;
    my_ylabel = ('Amplitude (dB (A))');
else
    my_ylabel = ('Amplitude (dB (L))');
end
figure(2),plot(freq,sensor_spectrum_dB,'b');grid
title(['Averaged FFT Spectrum  / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(freq(2)-freq(1)) ' Hz ']);
xlabel('Frequency (Hz)');ylabel(my_ylabel);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% display 3 : time / frequency analysis : spectrogram demo
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[sg,fsg,tsg] = specgram(signal,NFFT,Fs,hanning(NFFT),floor(NFFT*OVERLAP));  
% FFT normalisation and conversion amplitude from linear to dB (peak)
sg_dBpeak = 20*log10(abs(sg))+20*log10(2/length(fsg));     % NB : X=fft(x.*hanning(N))*4/N; % hanning only
% apply A weigthing if needed
if option_w == 1
    pondA_dB = pondA_function(fsg);
    sg_dBpeak = sg_dBpeak+(pondA_dB*ones(1,size(sg_dBpeak,2)));
    my_title = ('Spectrogram (dB (A))');
else
    my_title = ('Spectrogram (dB (L))');
end
% saturation of the dB range : 
% saturation_dB = 60;  % dB range scale (means , the lowest displayed level is XX dB below the max level)
min_disp_dB = round(max(max(sg_dBpeak))) - spectrogram_dB_scale;
sg_dBpeak(sg_dBpeak<min_disp_dB) = min_disp_dB;
% plots spectrogram
figure(3);
imagesc(tsg,fsg,sg_dBpeak);colormap('jet');
axis('xy');colorbar('vert');grid
title([my_title ' / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(fsg(2)-fsg(1)) ' Hz ']);
xlabel('Time (s)');ylabel('Frequency (Hz)');
function pondA_dB = pondA_function(f)
	% dB (A) weighting curve
	n = ((12200^2*f.^4)./((f.^2+20.6^2).*(f.^2+12200^2).*sqrt(f.^2+107.7^2).*sqrt(f.^2+737.9^2)));
	r = ((12200^2*1000.^4)./((1000.^2+20.6^2).*(1000.^2+12200^2).*sqrt(1000.^2+107.7^2).*sqrt(1000.^2+737.9^2))) * ones(size(f));
	pondA = n./r;
	pondA_dB = 20*log10(pondA(:));
end
function  [freq_vector,fft_spectrum] = myfft_peak(signal, Fs, nfft, Overlap)
% FFT peak spectrum of signal  (example sinus amplitude 1   = 0 dB after fft).
% Linear averaging
%   signal - input signal, 
%   Fs - Sampling frequency (Hz).
%   nfft - FFT window size
%   Overlap - buffer overlap % (between 0 and 0.95)
signal = signal(:);
samples = length(signal);
% fill signal with zeros if its length is lower than nfft
if samples<nfft
    s_tmp = zeros(nfft,1);
    s_tmp((1:samples)) = signal;
    signal = s_tmp;
    samples = nfft;
end
% window : hanning
window = hanning(nfft);
window = window(:);
%    compute fft with overlap 
 offset = fix((1-Overlap)*nfft);
 spectnum = 1+ fix((samples-nfft)/offset); % Number of windows
%     % for info is equivalent to : 
%     noverlap = Overlap*nfft;
%     spectnum = fix((samples-noverlap)/(nfft-noverlap));	% Number of windows
    % main loop
    fft_spectrum = 0;
    for i=1:spectnum
        start = (i-1)*offset;
        sw = signal((1+start):(start+nfft)).*window;
        fft_spectrum = fft_spectrum + (abs(fft(sw))*4/nfft);     % X=fft(x.*hanning(N))*4/N; % hanning only 
    end
    fft_spectrum = fft_spectrum/spectnum; % to do linear averaging scaling
% one sidded fft spectrum  % Select first half 
    if rem(nfft,2)    % nfft odd
        select = (1:(nfft+1)/2)';
    else
        select = (1:nfft/2+1)';
    end
fft_spectrum = fft_spectrum(select);
freq_vector = (select - 1)*Fs/nfft;
end

MATLAB: Usinf fft shift i need only to display the positive part of the spectrum

Try this:

fs=1e5; %Sampling Frequency
t=0:1/(fs):0.001; %Time at which signal displays
%First,generating a signal
f=1e4;%Signal is 10kHz
Fn = f/2;                                                               % Nyquist Freqency (Added)
amp_signal=1;
signal=amp_signal*(sin(2*pi*f*t));
L = length(signal);                                                     % Added
%Frequency Domain of transmitted signal
transmited_signal_fft=fft(signal);
f_axis = linspace(0, 1, fix(L/2)+1)*Fn;                                 % Changed

Iv = 1:length(f_axis);                                                  % Index Vector (Added)
% % transmited_signal_fft=fftshift(transmited_signal_fft);              % (Commented-Out)
% % f_axis=linspace(-fs/2,fs/2,length(transmited_signal_fft));          % (Commented-Out)
%Plotting Signal
figure;
subplot(2,1,1);
plot(t,signal);
grid on;
title('Transmitted signal in time domain')
%Display transmitted signal in frequency domain
subplot(2,1,2);
plot(f_axis, 2*abs(transmited_signal_fft(Iv))/L);                       % Changed
grid on;
title('Transmitted signal in frequency domain')

See this version of the documentation for details: fft (link).

Best Answer

Related Solutions

MATLAB: Time domain and frequency domain signal generation and plot

MATLAB: Usinf fft shift i need only to display the positive part of the spectrum

Related Question