Tag: signal processing

DSP Lab – Week 1

~ Kishore Athrasseri ~ 1 Comment

Constructing the Complex Plane

Suppose we have a sampled signal defined by the sequence h(n), n=0,1,2,...,N-1

Its Z- transform is given by H(z) = sum_{n=0}^{N-1} h(n)z^{-n} .

It maps the original sequence into a new domain, which is the complex plane z=e^{sT} where s=sigma+jomega is the parameter in the Laplace domain and T is the sampling period.

The jomega axis in the s-plane maps onto the unit circle with centre at the origin in the z-plane. So the value of H(z) at different points on the unit circle actually gives the contribution of the frequency component given by angle z, in the original signal.

This, in effect, gives the Discrete Fourier Transform of the sequence. Consider the following example:

%original sequence
h = [1,2,3,4];

%number of chosen points on the unit circle
N = 64;

%define the chosen points
z = complex(cos(2*pi/N*(0:N-1)),sin(2*pi/N*(0:N-1)));

%evaluate H(z) at each point
for i = 1:N
H(i) = 1+2*z(i)^-1+3*z(i)^-2+4*z(i)^-3;
end

%plot the unit circle
plot(z)

%plot the value of H(z) along the unit circle
figure
plot(abs(H))

%plot the N-point DFT of h(n)
figure
plot(abs(fft(h,64)))

This example computes the value of H(z) at 64 uniformly spaced points on the unit circle and compares it with the 64 point DFT. We can see that both (fig. b & c) are identical.

lab1_a
unit circle
lab1_b
value of H(z)
abs(fft(h))
abs(fft(h))

FFT with Phoenix

~ Kishore Athrasseri ~ Leave a comment

A companion analog box for Phoenix, is being developed here at IAUC, with which you can study AM, FM etc. It’s still in the process of development, but I had observed AM and FM signals generated by it, and studied their frequency spectra using the fft() in SciPy. We felt that it would be convenient to add these functions to the Python library of Phoenix, as given below. You can use these functions by adding this code segment to the file /usr/lib/python2.4/site-packages/phm.py

But since this bit of code has not really been tested, it’s better to avoid mixing these functions with the original phm.py . Instead, you can save this file to your home directory, and import it along with phm. It can be used just like phm, and you can initialize an object by typing “s=phsig.phsig()”.

    def fft(self, data):
        """
        Returns the Fourier transform of the signal represented by the samples
        in 'data' which is obtained by a read_block() call.
        Usage example:
        x = p.read_block(200,10,1)
        y = p.fft(x) 
        """
        from scipy import fft       
        np = len(data)
        delay = data[1][0] - data[0][0]

        v = []
        for i in range(np):
            v.append(data[i][1])

        ft = fft(v)

        #corrections in the frequency axis
        fmax = 1/(delay*1.0e-6)
        incf = fmax/np
        freq = []
        for i in range(np/2):
            freq.append(i*incf - fmax/2)
        for i in range(np/2,np):
            freq.append((i-np/2)*incf)

        ft_freq = []

        for i in range(np/2):
            x = [freq[i], abs(ft[i+np/2])/np]
            ft_freq.append(x)    

        for i in range(np/2,np):
            x = [freq[i], abs(ft[i-np/2])/np]
            ft_freq.append(x)
        return ft_freq                

    def plot_fft(self, data):              
        """
        Plots the Fourier transform of the signal represented by the samples in
        'data'. Calls self.fft() for obtaining the Fourier Transform
        Usage example:
        x = p.read_block(200,10,1)
        p.plot_fft(x)
        """
        ft_freq = self.fft(data)
        np = len(data)
        delay = data[1][0] - data[0][0]
        fmax = 1/(delay*1.0e-6)

        y = 0.0
        for i in range(np):
            if ft_freq[i][1] > y:
                y = ft_freq[i][1]

        if self.root == None:
            self.window(400,300,None)
        self.remove_lines()
        self.set_scale(-fmax/2, 0, fmax/2, y*1.1)
        self.line(ft_freq)

    def freq_comp(self, data):
        """
        Returns and displays the frequency components with the greatest
        spectral density. Calls self.fft() for obtaining the Fourier transform
        Usage example:
        x = p.read_block(200,10,1)
        p.freq_comp(x)    #Only prints the components
        y = p.freq_comp(x) #Prints and stores the components in a variable
                           #as an array of [strength, component]     
        """
        ft_freq = self.fft(data)
        np = len(data)
        delay = data[1][0] - data[0][0]

        peaks = []
        for n in range(1,np-1):
            a = ft_freq[n-1][1]
            b = ft_freq[n][1]
            c = ft_freq[n+1][1]           
            if (b>50) & (b>a) & (b>c):
                peaks.append([ft_freq[n][1],ft_freq[n][0]])

        peaks.sort()
        peaks.reverse()
        print 'Dominant frequency components are:'
        for i in range(len(peaks)):
            print '%6.3f Hz, %5.3f mV'%(peaks[i][1],peaks[i][0])
        return peaks