Hi all,

I am trying to model spectrum of a single rectangular pulse in matlab. 

I am using sinc formula versus fft and got some two issues that I can't understand.. 

issue 1: For pulse duration (T) > 1 there is small difference between fft result and sinc formula that gets better as I increase T.

issue 2: for T = 1 the sinc sees it correctly as unitary rectangular pulse but fft sees it as impulse since I set only one sample to 1. 

for T = 0 the sinc formula sees it as impulse understandably. fft would be dead

In both methods I use same discrete frequency grid.

Is this expected? Any contribution appreciated.

Here is my matlab code, (try changing T value)

___________________________________

n = 1024;                       %frequency grid resolution

T = 10;                           % pulse duration

f = freqspace(n,'whole'); %full circle frequency grid 

f = f/2 - .5;                      %for sinc must be -/+ rather than 0~1

%sinc method:

X = sinc(T*f);                  %same as sin(pi*T*f)/(pi*T*f)       

plot(f,abs(X),'-');              %amplitude cannot be negative

hold;

%fft methd

y = fft(ones(1,T),n)/T;

y = fftshift(y);

plot(f,abs(y),'g--');

It looks like you've already been pointed in the right direction, but I'll add that this confusion is often exacerbated by the casual use of the term "sinc" to describe anything reasonably sinc-ish shaped.   This is motivated by the infinitely large number of possible sinc-ish shapes that are encountered in filter design and windowing and antenna design and a number of other places.

It is practical and common to use the term that way and it often is, and, as happens often with all kinds of things like this, if you know to ask for clarification if you think "sinc" means "sinc-ish" and not definitively a sinc, it is good to ask clarifying questions.

And when you're writing (or talking) about it, it is often a judgement call whether to use "sinc" or some other cumbersome description that is overly specific when it doesn't need to be.

Anyway, I just said all that to say that this is sometimes, or often, the source of the confusion when things like this come up.  This is sometimes how people learn the difference between a sinc and a Dirichlet kernel and the other infinite number of similar shapes that might be encountered. 

If I’m understanding correctly, that make intuitive sense. In continuous time, the Fourier transform of a rectangular pulse is a sinc function. If the sinc function hasn’t died out before the Nyquist frequency, then there would be aliasing. You can use the transfer function of a recursive moving average filter to figure out what the actual FFT bins are.

Hi dszabo,

Thanks. It does make sense regarding aliasing and may explain the difference but that is if the rectangular pulse is in analogue world. I thought the fft being in digital domain will only check frequencies between -.5 and +.5 and any higher frequency cannot exist in my digital model or can it?

That is true, but the sinc function would only be correct to use in continuous time, but we can still use it as a model for why we get a different result in discrete time. If we want to establish a numerical solution in the discrete time domain, the easiest way I can think of is to calculate the z-transform of a rectangular pulse in discrete time. I recommend using the recursive moving average because it is easier to work with.  The engineers and scientists guide to dsp provides a solution in its chapter on moving average filters and it should be something like: sin(Nk)/sin(k). Note that for N=1, we get unity, and for large N we get something like a periodic sinc function.

My apologies, a moving average filter is a rectangular pulse, so the mathematics of deriving the frequency response of the moving average filter are identical to finding the Fourier transform of a rectangular pulse

Many Thanks I now got it right. It seemed there are different versions of rectangular pulse equations. The one that works is sin(T*pi*f)./(T*sin(pi*f));

But strangely it is not exactly sin(x)/x

here is the corrected code:

n = 1024;                       %frequency grid resolution

T = 10;                           % pulse duration

f = freqspace(n,'whole'); %full circle frequency grid 

f = f/2 - .5;                      %for sinc must be -/+ rather than 0~1

%sinc method:

X = sin(T*pi*f)./(T*sin(pi*f));   %this is the change  

plot(f,abs(X),'-');              %amplitude cannot be negative

hold;

%fft methd

y = fft(ones(1,T),n)/T;

y = fftshift(y);

plot(f,abs(y),'g--');

Glad to hear you got it sorted!  Remember that sin(x) is approximately x when x is small.  That is why it mostly looks like a sinc function for small values and not at high values

That would be equation (27) in my article I referenced (adjusted to T).

Ced

Hi again,

I am a bit clueless as what has averaging filter got to do here. Do you mean applying it to the pulse before sinc frmula or applying it to the pulse before fft or do you have other ideas. After all I want to analyse the pulse and not filter it.

Thanks

Actually, on re-reading, it is probably merited to point out that the frequency content repeats infinitely.  This is why the z-domain is represented in polar coordinates.  Where z=r*e^jw, we would get the same value for w=x that we would for w=x+2pi. It’s not that the higher frequency content doesn’t exist, it’s that it’s redundant

The sinc function is for the continuous case and is only an approximation in the discrete case.

It is fairly straightforward to derive the DFT bin values of a rectangular pulse function.  See the math in my blog article: "DFT Bin Value Formulas for Pure Complex Tones" which is on this site at:

dsprelated.com/showarticle/1038.php

In equation (9) set alpha and phi equal to zero, and modify the range of the summation to be over your pulse region.  The do a change in variable to convert the summation from being 0 to T-1 (where T is your pulse length).  The rest will follow naturally.

If you are stumped by this, I'll do the math and post it here.  The bottom line is the proper function is not the sinc function, but what is known as either the alias sinc function, or Dirichlet kernel.

Hope this helps,

Ced

Hi Ced,

Thanks for the response. What about DAC sinc envelope due to zero hold. Does that follow a continuous sinc formula which I think is sin(pi*f/Fs)/(pi*f/Fs). 

Yes.  The formal method that I learned for modeling mixed continuous- and discrete-time systems is to model the discrete-time signals as trains of impulses with weights equal to the actual signal, and model a DAC as a continuous-time filter that's rectangular in time -- which means it has a perfect sinc function in the frequency domain.

Hi kaz,

That is a different issue, one I am not that knowledgeable in, but I believe that is correct.

Ced

Chapter 3 of the Engineers and Scientists Guide to DSP!  A fantastic (and free) publication. You are correct

You are misapplying the math.  You are comparing apples to oranges.  You are -- well, you are making words fail me.

The sync function is the Fourier transform of a continuous time rectangular pulse happening in infinite time.  The FFT is a DFT, which is not just a Fourier transform of discrete-time data, but it is the Fourier transform of discrete-time data that comes in finite-sized chunks.  In Fourier-land, you couldn't be further away and still be dealing with single-dimensioned real data.

The DFT (and thus the FFT) is exact, but when you sample a rectangular pulse you are making an approximation (you are, in fact, generating a huge amount of energy in aliases).  When you limit the time domain to a finite interval, you are making another approximation.  It should be no surprise at all that you get different answers from the two different cases.

yes that was true in my first code but corrected later in subsequent post. The corrected code is given in that post.

Here is the problem...

The continuous rectangle spanning the interval -NT to +NT has a continuous transform often called the sinc of the form h(f)= 1/(2*NT) *sin(2*pi*f*2NT/2)/(2*pi*f*2NT/2)

the sampled data version of the rectangle containing 2N+1 samples T seconds apart has a continuous periodic transform called the Dirichlet Kernel (The periodic extension of the sinc function with spectral spacing matching the sample rate fs =1/T.

Its form is sin((2N+1)*theta/2)/sin(theta/2). For small angles the Dirichlet kernal approximates the sinc. The sinc spectrum falls like 1/f while the periodic extension crosses the half sample rate with zero slope due to the interaction of the periodic replicas.

The zeros of the sinc reside on the infinite line while the zeros of the Dirichlet reside on the circle.

If N is an even integer, say 10 for instance, when we project the alternating signs of cos(pi n) on the sampled rectange we would find the sum = +1 which says the Dirichlet crosses the half sample rate (at pi) with sidelobe value = +1; with the peak value at theta = 0 is 2N+1 or 21

If N is an odd integer, say 11 for instance, when we project the alternating signs of cos(pi n) on the sampled rectange we would find the sum = -1 which says the Dirichlet crosses the half sample rate (at pi) with sidelobe value = -1;   with the peak value at theta = 0 is 2N+1 or 23 

If we actually sample the continuous fectangle we have the problem of sampling the pulse at its discontinuities at -NT and +NT. The sample set would then contain 2N-1 ones bracketed by samples 0.5... of the form [0.5 1 1 1 .... 1 0.5]; this is the same as a rectangle sample set of 2N-1 samples convolved wih [0.5 0.5] which forms a windowed Dirichlet kernel...  with zeros at the crossover frequency theta=pi. All three options are illustrated in the attached matlab file     

f_h1a=zeros(1,211);f_h1a(107+(-10:10))=ones(1,21); 
f_h1=fftshift(f_h1a);

figure(1);subplot(2,1,1);
plot((-0.5-1/422:1/211:0.5-2/422)*211,f_h1a,'linewidth',2)
grid on;axis([-15 15 -0.1 1.2])
title('Spectrum, -10:10, ones(1,21) in 211 Sample Array')
xlabel('Frequency Index')
ylabel('Amplitude')

subplot(2,1,2)
plot(-105:105,fftshift(real(fft(f_h1))),'linewidth',2);hold on
plot(-100:10:-10,zeros(1,10),'ro','linewidth',2)
plot(10:10:100,zeros(1,10),'ro','linewidth',2)
plot(0, 21, 'ro','linewidth',2);hold off;grid on
axis([-106 106 -7 23]);
title('Time Response, Dirichlet Kernal, sin(21*\theta/2)/sin(\theta/2)')

xlabel('Time Index');ylabel('Amplitude');
f_h2a=zeros(1,231);f_h2a(117+(-11:11))=ones(1,23);
f_h2=fftshift(f_h2a);

figure(2);subplot(2,1,1)
plot((-0.5-1/462:1/231:0.5-2/462)*231,f_h2a,'linewidth',2)
grid on;axis([-15 15 -0.1 1.2])
title('Spectrum, -11:11, ones(1,23) in 231 Sample Array')
xlabel('Frequency Index');ylabel('Amplitude')
subplot(2,1,2)
plot(-115:115,fftshift(real(fft(f_h2))),'linewidth',2);hold on
plot(-110:10:-10,zeros(1,11),'ro','linewidth',2)
plot(10:10:110,zeros(1,11),'ro','linewidth',2)
plot(0, 23, 'ro','linewidth',2);hold off;grid on;
axis([-116 116 -7 25])
title('Time Response, Dirichlet Kernal, sin(23*\theta/2)/sin(\theta/2)');xlabel('Time Index');
ylabel('Amplitude');f_h3a=zeros(1,211);
f_h3a(107+(-10:10))=[0.5 ones(1,19) 0.5];f_h3=fftshift(f_h3a);
figure(3);subplot(2,1,1)
plot((-0.5-1/422:1/211:0.5-2/422)*211,f_h3a,'linewidth',2)
grid on;axis([-15 15 -0.1 1.2])
title('Spectrum, [-10:10], [0.5 ones(1,19) 0.5] in 211 Sample Array');xlabel('Frequency Index');ylabel('Amplitude')
subplot(2,1,2);
plot(-105:105,fftshift(real(fft(f_h3))),'linewidth',2);hold on
plot(-105:10.5:-10,zeros(1,10),'ro','linewidth',2)
plot(10:10.5:105,zeros(1,10),'ro','linewidth',2)
plot(0, 20, 'ro','linewidth',2);hold off;grid on;
axis([-116 116 -7 25])
title('Time Response, Dirichlet Kernal, sin(19*\theta/2)/sin(\theta/2)* cos(\theta/2)')
xlabel('Time Index');ylabel('Amplitude')

I couldn't attach matlab file so here it is in text

f_h1a=zeros(1,211);

f_h1a(107+(-10:10))=ones(1,21);

f_h1=fftshift(f_h1a);

figure(1)

subplot(2,1,1)

plot((-0.5-1/422:1/211:0.5-2/422)*211,f_h1a,'linewidth',2)

grid on

axis([-15 15 -0.1 1.2])

title('Spectrum, -10:10, ones(1,21) in 211 Sample Array')

xlabel('Frequency Index')

ylabel('Amplitude')

subplot(2,1,2)

plot(-105:105,fftshift(real(fft(f_h1))),'linewidth',2)

hold on

plot(-100:10:-10,zeros(1,10),'ro','linewidth',2)

plot(10:10:100,zeros(1,10),'ro','linewidth',2)

plot(0, 21, 'ro','linewidth',2)

hold off

grid on

axis([-106 106 -7 23])

title('Time Response, Dirichlet Kernal, sin(21*\theta/2)/sin(\theta/2)')

xlabel('Time Index')

ylabel('Amplitude')

f_h2a=zeros(1,231);

f_h2a(117+(-11:11))=ones(1,23);

f_h2=fftshift(f_h2a);

figure(2)

subplot(2,1,1)

plot((-0.5-1/462:1/231:0.5-2/462)*231,f_h2a,'linewidth',2)

grid on

axis([-15 15 -0.1 1.2])

title('Spectrum, -11:11, ones(1,23) in 231 Sample Array')

xlabel('Frequency Index')

ylabel('Amplitude')

subplot(2,1,2)

plot(-115:115,fftshift(real(fft(f_h2))),'linewidth',2)

hold on

plot(-110:10:-10,zeros(1,11),'ro','linewidth',2)

plot(10:10:110,zeros(1,11),'ro','linewidth',2)

plot(0, 23, 'ro','linewidth',2)

hold off

grid on

axis([-116 116 -7 25])

title('Time Response, Dirichlet Kernal, sin(23*\theta/2)/sin(\theta/2)')

xlabel('Time Index')

ylabel('Amplitude')

f_h3a=zeros(1,211);

f_h3a(107+(-10:10))=[0.5 ones(1,19) 0.5];

f_h3=fftshift(f_h3a);

figure(3)

subplot(2,1,1)

plot((-0.5-1/422:1/211:0.5-2/422)*211,f_h3a,'linewidth',2)

grid on

axis([-15 15 -0.1 1.2])

title('Spectrum, [-10:10], [0.5 ones(1,19) 0.5] in 211 Sample Array')

xlabel('Frequency Index')

ylabel('Amplitude')

subplot(2,1,2)

plot(-105:105,fftshift(real(fft(f_h3))),'linewidth',2)

hold on

plot(-105:10.5:-10,zeros(1,10),'ro','linewidth',2)

plot(10:10.5:105,zeros(1,10),'ro','linewidth',2)

plot(0, 20, 'ro','linewidth',2)

hold off

grid on

axis([-116 116 -7 25])

title('Time Response, Dirichlet Kernal, sin(19*\theta/2)/sin(\theta/2)* cos(\theta/2)')

xlabel('Time Index')

ylabel('Amplitude')