How to Highlight the Source Code in Blogger:

I had been through a lot for searching the solution to highlight the source code in my Blogger.

This was also the reason I give up on the use of Blogger.

Then I move to here, which has a fine support of highlighting source code;

Here is the step to follow:

200 things to do with Linux

Copyright @ 200th issue of Linux Journal.
1. Actually work instead of waiting for reboots.—Tim Chase
2. Add extra monitors.—LJ Staff
3. Analyse water level and precipitation data.—Keith Nunn
4. Analysis of remote sensing imagery.—Micha Silver
5. Antagonize Windows users.—John Abbott
6. Anything I need, since 1994.—Manuel Trujillo
7. As the basis for FOSS conferences.—moose
8. Audio chat.—LJ Staff
9. Automate tasks with bash.—Dusty Roberson
10. Avoid using Microsoft Windows!—Simon Quantrill, Chris Szilagyi

Network Programming Project #03

Description:

The purpose of this project will be to simulate the performance of a first-in-first-out (FIFO)
queue with fixed sized packets and Markov arrivals to the queue. The simulation will be
performed using only one input source to the queue.

The following parameters will be used throughout the entire project:
c = link speed = 10 Mbps
p = packet size = 4000 bits

Simulation of a FIFO Queue with a Single Source:

Create your own, newly created simulation program for simulation of the queue.
You should write a program that does the low-level management of an event list as described in
the lecture.

Your purpose will be to simulate the performance of the queue for various packet arrival rates.
Within one simulation, you must simulate many packets and find the average response time for
those packets. Then you should perform multiple iterations of the simulation so as to compute
the average response times for packets through the M/D/1 queue across several simulation runs.

Use the following parameters:
λ = arrival rates = 1000 to 2000 in steps of 100
10 simulation runs per arrival rate to compute the average response times
The simulation time should be long enough to simulate 1000 packet arrivals.

To generate the exponential random number needed to create the interarrival times for the
packets, you can use your programming environment’s random number generator. Use the
function that generates a random number X uniformly between 0 and 1. Then use the following
formula to create an exponentially distributed random number Y, where ln(X) is the natural
logarithm function.

Problem 01:

A single plot that includes the response time for each simulation at each arrival rate,
the average response time at each arrival rate, and the theoretical response time as
given by the following equation.

#!/usr/bin/env  python
# Copyright (C) 2010-11-13 by Antonio081014
 
import random
import numpy
import matplotlib.pyplot
 
LinkRate = 10000000 #kbps;
pkgSize = 4000      #bits;
simTimes = 10       #Number of simulations;
longda = numpy.arange(1000, 2001, 100)
pkgNum = 2000
servTime = 1. * pkgSize / LinkRate
M = servTime
def getArrivalTime(lnda):
    r = random.random()
    return -numpy.log(r) / lnda;
 
def performance(rate):
        startTime = 0.
    deptTime = servTime
    totalTime = servTime
    for i in range(1, pkgNum):
        startTime += getArrivalTime(rate)
        if deptTime - startTime <= 0.:
            deptTime = startTime + servTime
            totalTime += servTime
        else:
            totalTime += deptTime + servTime - startTime;
            deptTime += servTime
    return 1. * totalTime / pkgNum
 
def simulate():
    ret = numpy.zeros((simTimes, longda.shape[0]), 'float')
    for count in range(simTimes):
        for r in range(longda.shape[0]):
            rate = longda[r]
            ret[count][r] = performance(rate);
    return (ret)*1000.
 
def getTheory(lnda):
    ld = numpy.zeros(lnda.shape[0], 'float')
    for i in range(lnda.shape[0]):
        d = lnda[i]
        temp = (2.*M - d*(M**2))/(2.*(1.-d*M));
        ld[i] = temp
    return ld*1000.
 
def plotSimulate(record):
    m = record.shape[0]
    n = record.shape[1]
    sum = numpy.sum(record, axis=0);
    temp = numpy.arange(n)*100 + 1000
    matplotlib.pyplot.plot(temp, sum/m, 'g', label='Simulated Average Response Time');
    matplotlib.pyplot.plot(temp, getTheory(longda), 'r', label='Theoretical Average Response Time');
    for i in range(n):
        for j in range(m):
            matplotlib.pyplot.plot(i*100+1000, record[j][i], 'r*');
    matplotlib.pyplot.legend();
    matplotlib.pyplot.xlabel('Arrival Rate (packages per second)');
    matplotlib.pyplot.ylabel('Response Time (ms)');
    matplotlib.pyplot.title('Simulated Response Time for an M/D/1 Queue');
 
if __name__ == '__main__':
    record = simulate();
#   print record
    matplotlib.pyplot.ion()
    matplotlib.pyplot.figure()
    plotSimulate(record);
    matplotlib.pyplot.savefig('SimulationT.png')

Problem 02:

Two plots of actual queue fill with respect to time for one simulation. Show the queue
fill at the time of each packet arrival before the packet is entered into the queue or
serviced. The first plot should be for an arrival rate of 2000 packets per second, and
the second plot should be for an arrival rate of 3000 packets per second.

The result figure:

#!/usr/bin/env  python
# Copyright (C) 2010-11-13 by Antonio081014
# Single arrival rate specified
 
import random
import numpy
import matplotlib.pyplot
 
LinkRate = 10000000 #kbps;
pkgSize = 4000      #bits;
simTimes = 10       #Number of simulations;
longda = numpy.array([2000, 3000])
pkgNum = 1000
servTime = 1. * pkgSize / LinkRate
M = servTime
def getArrivalTime(lnda):
    r = random.random()
    return -numpy.log(r) / lnda;
 
def performance(rate):
    start = numpy.zeros((pkgNum), 'float');
    dept = numpy.zeros((pkgNum), 'float');
        startTime = 0.
    deptTime = servTime
    totalTime = servTime
    start[0] = startTime
    dept[0] = deptTime
    for i in range(1, pkgNum):
        startTime += getArrivalTime(rate)
        if deptTime - startTime <= 0.:
            deptTime = startTime + servTime
            totalTime += servTime
        else:
            totalTime += deptTime + servTime - startTime;
            deptTime += servTime
        start[i] = startTime
        dept[i] = deptTime
    return start*1000., dept*1000.
 
def simCount(rate):
    start,dept = performance(rate)
    count = numpy.zeros((pkgNum));
    for i in range(pkgNum):
        for j in dept[:i]:
            if j > start[i]:
                count[i]+=1
    return count
     
 
def plotSimulate(rate):
    record = simCount(rate);
    m = record.shape[0]
    temp = numpy.arange(m)
    matplotlib.pyplot.plot(temp, record, 'r*', label='The situation of filled queue before next package\'s arrival');
    matplotlib.pyplot.legend();
    matplotlib.pyplot.xlabel('Package Number');
    matplotlib.pyplot.ylabel('Number of package in the queue when its comming');
    matplotlib.pyplot.title('The sistuation of filled queue before next package\'s arrival when arrival rate is ' + str(rate));
 
if __name__ == '__main__':
    matplotlib.pyplot.ion()
    matplotlib.pyplot.figure()
    plotSimulate(longda[0]);
    matplotlib.pyplot.savefig('Rate01.png')
    matplotlib.pyplot.figure()
    plotSimulate(longda[1]);
    matplotlib.pyplot.savefig('Rate02.png')

The result figures:

Matlab/Voice Session.

The Content For Matlab/Voice Session.
GuessNumber.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function GuessNumber
    disp('Generating a random number from 0 to 10');
    disp('You have 7 times to guess this number');
    disp('I will tell you if your guessed number is greater of less then the right number');
    number = ceil(rand(1,1)*10);
    count = 0;
    while(count < 5)
        temp = input('What Integer You Want to Try: ');
        if temp == number
            if count <= 1
                disp('Genieus, You got it right');
            else
                disp(['You are so smart, It is exactly ', num2str(temp)]);
            end
            break;
        else
            if temp > number
                disp(['Try some numbers smaller than ', num2str(temp)]);
            else
                disp(['Try some numbers larger than ', num2str(temp)]);
            end
            if count <=2
                disp('You still have chance to win, Good Luck');
            else
                disp([num2str(5-count), ' times left, Take care what you choose']);
            end
        end
        count = count + 1;
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
LeapLeap.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function LeapLeap
%     year = 2000;
    year = input('Type the year you want to check, e.g. 2000 : ');
    if (mod(year, 4)==0 && mod(year, 100)~=0) || mod(year, 400)==0
        disp([num2str(year), ' is a leap year']);
    else
        disp([num2str(year), ' is NOT a leap year']);
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Voice Part.
1st Links.
2nd Links.
3rd Links.

Speech Processing, Proj06

//Purdue Cal - ECE 595C.
//Spring 2010.
//Project 06.
//Copyright @ 2010 antonio081014 ;
//All codes are in Matlab.

Report Link:

Task:
1. Speech enhancement by using spectral subtraction;




%% Proj06;
clear all;
clc;
close all;

fd = fopen('chj.16', 'rb', 'b');
[spch, K] = fread(fd, 'int16');
fclose(fd);

% Remove the explosive in the background;
spch(28080:28100) = 0;
sspch = (spch - mean(spch));
y = sspch /max(sspch);

fs = 8000;
FL = 20 * fs / 1000; % 20ms;
NF = 1024;
frame = zeros(4, FL);
frame(1,:) = y(2241:2400);
frame(2,:) = y(6001:6160);
frame(3,:) = y(14912:15071);
frame(4,:) = y(33226:33385);

STFT = zeros(4, NF); % short-time spectrum;
STauto = zeros(4, 2 * FL -1); % short-time autocorrelation;
STPow = zeros(4, NF); % short-time
for i=1:4
    STFT(i, :) = fft(frame(i,:).*hamming(FL)', NF);
    STauto(i, :) = xcorr(frame(i,:));
    STPow(i, :) = abs(fft(frame(i,:).*hamming(FL)', NF)).^2; %Using hamming window to reduce the error.
%     STPow(i, :) = abs(fft(frame(i,:), NF)).^2;
end

AvergSTPow = sum(STPow, 1) / 4;

Newspch = zeros(1, length(y));
Fstart = 1;
Fend = Fstart + FL - 1;
overlap = FL / 2;
Fcnt = 0;
while(Fend <= length(y))
  
    data = y(Fstart:Fend).*hamming(FL);
    FFTdata = fft(data, NF);
  
    Re = max(abs((FFTdata)').^2 - AvergSTPow, 0);  %Get the Power of the spectral differences;
    phase = angle((FFTdata)'); % Keep the phase of the original speech;
    Newdata = sqrt(Re) .* (cos(phase) + i * sin(phase));
%    Newdata = sqrt(Re) .* exp(i * phase); % This way doesn't work. I don't know why.
    tmp = real(ifft(Newdata));
    Newspch(1, Fstart: Fend) = Newspch(1, Fstart: Fend) + tmp(1, 1:160);
%  don't use specified length when using ifft.
%     Newspch(1, Fstart: Fend) = Newspch(1, Fstart: Fend) + real(ifft(Newdata, FL)); % Reconstruct the speech with 50% overlap.
    Fstart = Fstart + overlap;
    Fend = Fstart + FL - 1;
    Fcnt = Fcnt + 1;
end

% soundsc(y, fs);
% soundsc(Newspch, fs);
figure; hold on;
plot(y, 'r'); plot(Newspch, 'g');

Speech Processing, Proj05

//Purdue Cal - ECE 595C.
//Spring 2010.
//Project 05.
//Copyright @ 2010 antonio081014  ;
//All codes are in Matlab.


Report Link:

Task:
1. Getting the Mel-Scale Frequencies Cepstral Coefficient.
    1.1 Using overlapping Triangular window.
    1.2 Calculate the energy for each filter band. (16 in total).
    1.3 Set all values (except energy values on these 16 frequencies) to be zeros.
    1.4 Do Discrete Cosine Transform.
   Plot the first feature through all the frames for one utterance. (trajectories).

2. Dynamic Time Warping with simple Euclidean and cosine distance.
    Here, I didn't add any constraints in DTW.

3. Getting the features for equal spaced frequencies' magnitude.

    Use equal spaced frequencies is sure not better than Mel Scale frequencies. This is why we use MFCCs instead of using these. MFCCs can better approximate the human auditory system.

Speech Processing, Proj04

//Purdue Cal - ECE 595C.
//Spring 2010.
//Project 04.
//Copyright @ 2010 antonio081014  ;
//All codes are in Matlab.


Report Link:


Task one:
1. Get the Real Cepstrum of the windowed speech.
2. Separate the Vocal Tract and the Glottal part by using a low time lifter.
3. Do DFT to the h[n] in quefrency.
4. Get the H(w) by using exp(reverse of the operation 'log'). 


Hint:
1. How to read pitch frequency: use two of the adjacent excitation peaks.
2. How to read F1: from the result of step 4 above, we can easily read the F1, F2...(Formants).
3. Pitch frequency has nothing relationship with H(w) of VT(Vocal Tract).


Task Two:
1. Using DFT to get the spectrum of H(w);
2. Using LPC to get the envelop of the spectrum of H(w).
3. Using Cepstrum to get the spectrum of H(w).
4. Compare the F1, F2,... and also pitch frequency with each other method.


Task Three:
1. Using a unvoiced frame(fricative) to do the cepstrum.
2. Using a voiced frame to do the cepstrum.
3. Check the difference of F1, F2 ... with these two frames.

Speech Processing, Proj02

//Purdue Cal - ECE 595C.
//Spring 2010.
//Project 02.
//Copyright @ 2010 antonio081014  antonio081014 ;
//All codes are in Matlab.


Report Link:

Task 1: Obtain estimated F0  by using the log harmonic product spectrum with K = 4 and K=10.

PS: The log product of the harmonic spectrum is the summation of the log harmonic spectrum.\

1. Down sampling with K=4 and K=10 to get the spectrum for each case.

2. Sum them up.

3. Find the highest peak in the final spectrum. That is the estimated F0.

Code:

fdate = fft(dat(:,15).*hamming(320), 1024);  % The 15th frame in the data.

tmp = log(abs(fdat3(1:512)));

figure;

subplot(4,1,1);

plot(tmp);

title(['The frequency plot of the frame #:', num2str(23)]);

xlabel('Freq index');

ylabel('Log-Magnitude');

hold on;

tmp1 = downsample(tmp, 4);

z = zeros(512 - length(tmp1),1);

tmp1 = [tmp1; z];

tmp2 = downsample(tmp, 10);

z = zeros(512 - length(tmp2),1);

tmp2 = [tmp2; z];

subplot(4,1,2);

plot(tmp1, 'r');

xlabel('Freq index, K=4');

ylabel('Log-Magnitude');

subplot(4,1,3);

plot(tmp2, 'g');

xlabel('Freq index, K=10');

ylabel('Log-Magnitude');

product = tmp + tmp1 + tmp2;

subplot(4,1,4);

plot(product, 'r');

xlabel('Freq index');

ylabel('Log-Magnitude');

[x, y] = ginput(3);

% x = [11.79; 57.60; 93.18;]; These are the index in frequency domain.

Task 2: Using overlap-add method, reconstruct the given utterance after filtering each frame using a bandpass filter from 50Hz-2000Hz.

1. Filter each frame.

2. Reconstruct the utterance by using overlap-add method.

Code:

% Filter:

recdat = zeros(320, 100);

for i=1:100                                    % There is 100 frames here.

    fmdate = dat(:,i);

    ffmdat = fft(fmdate.*hamming(320),1024); % Take Fourier Transform using 1024 points.

    LF = ceil(1024 / fs * 50);

    HF = floor(1024 / fs * 2000);

    ffmdat(1:LF,1) = 0;

    ffmdat(HF:512,1) = 0;

    ffmdat(512:end-HF,1) = 0;

    ffmdat(end-LF:end,1) = 0;

    recdat(:,i) = real(ifft(ffmdat,320));

    clear ffmdat;

end

% Reconstruct

update_sp = zeros(1,16000);

start = 1;

stop = start + 320 - 1;

for i=1:100

    update_sp(start:stop) = update_sp(start:stop) + recdat(:,i)';    % overlap-add;

    start = start + 160;

    stop = start + 320 - 1;

    if stop > 16000

        break

    end

end

figure;

plot(update_sp);

title('The filterd speech.');

% soundsc(update_sp);

fm = recdat(:,23);

ffm = fft(fm.*hamming(320), 1024);

figure;

plot(abs(ffm(1:512)));

title(['The filterd frame #:', num2str(23)]);

[xx, yy] = ginput(3);

% xx = 35.2534562211982;72.5806451612903;107.142857142857;];

%%

figure;

x = nn_spch;

window = 128;

noverlap = 64;

nfft = 128;

fs = 16000;

spectrogram(x,window,noverlap,nfft,fs, 'yaxis');

colormap(gray);

title('Power Spectral Density, Original');

xlabel('In Time Domain.');

ylabel('In Freq Domain.');

figure;

x = update_sp;

window = 128;

noverlap = 64;

nfft = 128;

fs = 16000;

spectrogram(x,window,noverlap,nfft,fs, 'yaxis');

colormap(gray);

title('Power Spectral Density, Filterd');

xlabel('In Time Domain.');
ylabel('In Freq Domain.');

Speech Processing, Proj03

//Purdue Cal - ECE 595C.
//Spring 2010.
//Project 03.
//Copyright @ 2010 antonio081014 ;
//All codes are in Matlab.


Report Link:


Task One: How to plot the frequency spectrum and the spectrum envelop.
1. use the Matlab function lpc to generalize the coefficients alphas for the windowed samples in each frame.
2. use the Matlab function freqz to plot the spectrum envelop.
3. use the Matlab fft to generalize the features in Frequency domain, and then plot it in half (symmetry).
Code:

    for i=1:3
        ftmp = fft(data(:,i), 2048);   % data is the windowed signal.
        [a, g] = lpc(data(:,i), 7);       % 7th order.
        figure;
        hold on;
        [h, w] = freqz([1], a, 1024, fs);  % plot the spectrum envelop.
        plot(abs(h) / max(abs(h)), 'r');
        plot(abs(ftmp(1:1024)) / max(abs(ftmp(1:1024))), 'g'); % plot the freq spectrum.
        title(['Frame: #', num2str(i)]);
        legend('The spectrum envelop','The freq spectrum');
   end

It shows the spectrum envelop in red, and the frequency feature in green.
















Task Two: Generate the LP model by using Toeplitz Method.
1. Get the autocorrelation of the original signals.
2. Generate the matrix by using Matlab function Toeplitz .
3. Get the coefficients alphas according to the normal equations in the matrix form.

Code:
    tmp_data = norm_spch(201:600);
    c = xcorr(tmp_data,4, 'coeff');    % It can generate 2 * N -1 poles there. 4 -> 7.
    T = toeplitz(c(2:end));  % The first is 1, which is not necessary.
    alpha = inv(T) * c(2:end);
    lpca = lpc(tmp_data, 7);
    % Verified through the following steps.
    test_sp = filter([0  -lpca(2:end)], 1, tmp_data);
    plot(test_sp, 'r');
    hold on
    plot(tmp_data, 'g');
    err_sig = tmp_data - test_sp;
    figure;
    plot(err_sig);
    title('The error signal generated by using Toeplitz Method');
    xlabel('Time');
    ylabel('Magnitude');

It plots the original signal (green) and the estimated signal (red) generated by using Toeplitz method.