Signal Processing
In this page, I would post a quick reference for Matlab and Octave. (Octave is a GNU program which is designed to provide a free tool that work like Matlab. I don't think it has 100% compatability between Octave and Matlab, but I noticed that most of basic commands are compatible. I would try to list those commands that can work both with Matlab and Octave). All the sample code listed
here,
I tried with Octave, not with Matlab.
There are huge number of functions which has not been explained here, but I would try to list those functions which are most commonly used in most of matlab sample script you can get. My purpose is to provide you the set of basic commands with examples so that you can at least read the most of sample script you can get from here and there (e.g, internet) without overwhelming you. If you
get
familiar with these minimal set of functionality, you would get some 'feeling' about the tool and then you would make sense out of the official document from Mathworks or GNU Octave which explains all the functions but not so many examples.
I haven't completed 'what I think is the minimum set' yet and I hope I can complete within a couple of weeks. Stay tuned !!!
Signal Processing
< Quantization >
Ex)
Input
|
t = [0:.1:2*pi];
sig = sin(t);
partition = [-1:.2:1];
codebook = [-1.2:.2:1];
[index,quants] = quantiz(sig,partition,codebook);
plot(t,sig,'r-',t,quants,'b-',t,quants,'bo');
axis([-.2 7 -1.2 1.2]);
legend('Original signal','Quantized signal');
|
Output
|
|
< Zero Stuffing >
Ex 1)
Input
|
sig = [1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1];
sig_zerostuff = zeros(1,length(sig) * 2);
sig_zerostuff(1:2:length(sig_zerostuff)) = sig;
subplot(2,1,1);
stem(sig,'MarkerFaceColor',[0 0 1]);
subplot(2,1,2);
stem(sig_zerostuff,'MarkerFaceColor',[1 0 0]);
|
Output
|
|
Ex 2)
Input
|
N = 6; % Upsampling Rate
sig = [1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1];
sig_zerostuff = zeros(1,length(sig) * N);
sig_zerostuff(1:N:length(sig_zerostuff)) = sig;
subplot(2,1,1);
stem(sig,'MarkerFaceColor',[0 0 1]);
subplot(2,1,2);
stem(sig_zerostuff,'MarkerFaceColor',[1 0 0]);
|
Output
|
|
< Sample and Hold >
Ex 1)
Input
|
N = 6; % Upsampling Rate
h_t = ones(1,N);
sig = [1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1];
sig_zerostuff = zeros(1,length(sig) * N);
sig_zerostuff(1:N:length(sig_zerostuff)) = sig;
% this generate the last graph. If you are not sure how this routine can change the second graph
% to the third graph, refer to 'Convolution' page first.
sig_sample_hold = conv(sig_zerostuff,h_t);
subplot(3,1,1);
stem(sig,'MarkerFaceColor',[0 0 1]);
subplot(3,1,2);
stem(sig_zerostuff,'MarkerFaceColor',[1 0 0]);
subplot(3,1,3);
stem(sig_sample_hold,'MarkerFaceColor',[1 0 0]);xlim([1 length(sig_zerostuff)]);
|
Output
|
|
< Sample and Hold and Moving Avg >
Ex 1)
Input
|
N = 6; % Upsampling Rate
h_t = ones(1,N);
h_t_avg = (1/N) .* ones(1,N);
sig = [1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1];
sig_zerostuff = zeros(1,length(sig) * N);
sig_zerostuff(1:N:length(sig_zerostuff)) = sig;
sig_sample_hold = conv(sig_zerostuff,h_t);
sig_sample_hold_avg = conv(sig_sample_hold,h_t_avg);
subplot(4,1,1);
stem(sig,'MarkerFaceColor',[0 0 1]);
subplot(4,1,2);
stem(sig_zerostuff,'MarkerFaceColor',[1 0 0]);
subplot(4,1,3);
stem(sig_sample_hold,'MarkerFaceColor',[1 0 0]);xlim([1 length(sig_zerostuff)]);
subplot(4,1,4);
stem(sig_sample_hold_avg,'MarkerFaceColor',[1 0 0]);xlim([1 length(sig_zerostuff)]);
|
Output
|
|
< Sample and Hold and Low Pass Filter >
Ex 1)
Input
|
N = 6; % Upsampling Rate
h_t = ones(1,N);
h_t_lp = sinc(-pi:2*pi/20:pi); % you can use any of low pass filter coefficient you want
% I just tried to create a simple one without using signal
processing
% package
sig = [1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1];
sig_zerostuff = zeros(1,length(sig) * N);
sig_zerostuff(1:N:length(sig_zerostuff)) = sig;
sig_sample_hold = conv(sig_zerostuff,h_t);
sig_sample_hold_lp = conv(sig_sample_hold,h_t_lp);
subplot(4,1,1);
stem(sig,'MarkerFaceColor',[0 0 1]);
subplot(4,1,2);
stem(sig_zerostuff,'MarkerFaceColor',[1 0 0]);
subplot(4,1,3);
stem(sig_sample_hold,'MarkerFaceColor',[1 0 0]);xlim([1 length(sig_zerostuff)]);
subplot(4,1,4);
stem(sig_sample_hold_lp,'MarkerFaceColor',[1 0 0]);xlim([1 length(sig_sample_hold_lp)]);
|
Output
|
|
Ex 2)
: This example shows how "Interpolation" or "Upsampling" can be implemented by using 'zero stuffing' and 'low pass filter'. If you throw away the initial a couple of dozen samples of the last plot, you would see the sample plot as shown in 'interpolation' example.
Input
|
N = 4; % Upsampling Rate
h_t = ones(1,N);
h_t_lp = sinc(-pi:2*pi/20:pi);
h_t_lp = h_t_lp ./ sum(h_t_lp); % you can use any of low pass filter coefficient you want
% I just tried to create a simple one without using signal
processing
% package
Fs=1000;
A=1.5;
B=1;
f1=50;
f2=100;
t=0:1/Fs:1;
ti= 0:1/(Fs*N):1;
sig=A*cos(2*pi*f1*t)+B*cos(2*pi*f2*t);
sig_zerostuff = zeros(1,length(sig) * N);
sig_zerostuff(1:N:length(sig_zerostuff)) = sig;
sig_sample_hold = conv(sig_zerostuff,h_t);
sig_sample_hold_lp = conv(sig_sample_hold,h_t_lp);
subplot(4,1,1);
stem(sig,'MarkerFaceColor',[0 0 1]);xlim([1 25]);
subplot(4,1,2);
stem(sig_zerostuff,'MarkerFaceColor',[1 0 0]);xlim([1 25*N]);
subplot(4,1,3);
stem(sig_sample_hold,'MarkerFaceColor',[1 0 0]);xlim([1 25*N]);
subplot(4,1,4);
stem(sig_sample_hold_lp,'MarkerFaceColor',[1 0 0]);xlim([1 25*N]);
|
Output
|
|
Ex 3)
: This is exactly same as Ex 2 except that it has additional graph at the last step showing some initial samples removed to show the meaningfully interpolated graph.
Input
|
N = 4; % Upsampling Rate
h_t = ones(1,N);
h_t_lp = sinc(-pi:2*pi/20:pi);
h_t_lp = h_t_lp ./ sum(h_t_lp); % you can use any of low pass filter coefficient you want
% I just tried to create a simple one without using signal
processing
% package
Fs=1000;
A=1.5;
B=1;
f1=50;
f2=100;
t=0:1/Fs:1;
ti= 0:1/(Fs*N):1;
sig=A*cos(2*pi*f1*t)+B*cos(2*pi*f2*t);
sig_zerostuff = zeros(1,length(sig) * N);
sig_zerostuff(1:N:length(sig_zerostuff)) = sig;
sig_sample_hold = conv(sig_zerostuff,h_t);
sig_sample_hold_lp = conv(sig_sample_hold,h_t_lp);
subplot(5,1,1);
stem(sig,'MarkerFaceColor',[0 0 1]);xlim([1 25]);
subplot(5,1,2);
stem(sig_zerostuff,'MarkerFaceColor',[1 0 0]);xlim([1 25*N]);
subplot(5,1,3);
stem(sig_sample_hold,'MarkerFaceColor',[1 0 0]);xlim([1 25*N]);
subplot(5,1,4);
stem(sig_sample_hold_lp,'MarkerFaceColor',[1 0 0]);xlim([1 25*N]);
subplot(5,1,5);
stem(sig_sample_hold_lp,'MarkerFaceColor',[1 0 0]);xlim([1+length(h_t_lp)/2 25*N+length(h_t_lp)/2
|
Output
|
|
< Interpolation >
Ex)
Input
|
Fs=1000;
A=1.5;
B=1;
f1=50;
f2=100;
t=0:1/Fs:1;
ifreqs=4;
ti= 0:1/(Fs*ifreqs):1;
x=A*cos(2*pi*f1*t)+B*cos(2*pi*f2*t);
y=interp(x,ifreqs); % interpolate signal by ifreqs
stem(t(1:25),x(1:25),'MarkerFaceColor',[0 0 1]);hold on; % plot original signal
plot(ti(1:100),y(1:100),'ro','MarkerFaceColor',[1 0 0]);;hold off;
legend('original data','interpolated data');
|
Output
|
|
<interpolation - complex number >
Ex)
Input
|
clear all;
t = 0:pi/4:6*pi;
sig=exp(j*t);
intSig = zeros(1,2*length(sig));
intSig(1:2:end) = sig;
f = [0.0 0.3 0.6 1.0];
a = [1.0 1.0 0.0 0.0];
b = firpm(15,f,a); % this is filter design funtion with Remez Algorithm
[h,w] = freqz(b,1,16);
intFilterSig = 2 .* conv(intSig,b);
% plot original signal
subplot(4,4,1);
stem(abs(sig),'MarkerFaceColor',[1 0 0],'MarkerSize',3);
title('(1)');xlim([1 length(sig)]);ylim([-1.5 1.5]);
subplot(4,4,2);
stem(real(sig),'MarkerFaceColor',[1 0 0],'MarkerSize',3);
title('(2)');xlim([1 length(sig)]);ylim([-1.5 1.5]);
subplot(4,4,3);
stem(imag(sig),'MarkerFaceColor',[1 0 0],'MarkerSize',3);
title('(3)');xlim([1 length(sig)]);ylim([-1.5 1.5]);
subplot(4,4,4);
stem(angle(sig),'MarkerFaceColor',[1 0 0],'MarkerSize',3);
title('(4)');xlim([1 length(sig)]);ylim([-pi pi]);
% plot signal with zeros interleaved
subplot(4,4,5);
stem(abs(intSig),'MarkerFaceColor',[1 0 0],'MarkerSize',3);
title('(5)');xlim([1 length(intSig)]); ylim([-1.5 1.5]);
subplot(4,4,6);
stem(real(intSig),'MarkerFaceColor',[1 0 0],'MarkerSize',3);
title('(6)');xlim([1 length(intSig)]); ylim([-1.5 1.5]);
subplot(4,4,7);
stem(imag(intSig),'MarkerFaceColor',[1 0 0],'MarkerSize',3);
title('(7)');xlim([1 length(intSig)]); ylim([-1.5 1.5]);
subplot(4,4,8);
stem(angle(intSig),'MarkerFaceColor',[1 0 0],'MarkerSize',3);
title('(8)');xlim([1 length(intSig)]); ylim([-pi pi]);
% plot low pass filter
subplot(4,4,9);
plot(f,a,w/pi,abs(h));title('(9)');
subplot(4,4,10);
plot(w/pi,real(h));title('(10)');
subplot(4,4,11);
plot(w/pi,imag(h));title('(11)');
subplot(4,4,12);
stem(b,'MarkerFaceColor',[1 0 0],'MarkerSize',3);xlim([1 length(b)]);
title('(12)');
% plot signal with zeros interleaved and filtered
subplot(4,4,13);
stem(abs(intFilterSig),'MarkerFaceColor',[1 0 0],'MarkerSize',3);
title('(13)');xlim([1 length(intFilterSig)]); ylim([-1.5 1.5]);
subplot(4,4,14);
stem(real(intFilterSig),'MarkerFaceColor',[1 0 0],'MarkerSize',3);
title('(14)');xlim([1 length(intFilterSig)]); ylim([-1.5 1.5]);
subplot(4,4,15);
stem(imag(intFilterSig),'MarkerFaceColor',[1 0 0],'MarkerSize',3);
title('(15)');xlim([1 length(intFilterSig)]);ylim([-1.5 1.5]);
subplot(4,4,16);
stem(angle(intFilterSig),'MarkerFaceColor',[1 0 0],'MarkerSize',3);
title('(16)');xlim([1 length(intFilterSig)]); ylim([-pi pi]);
|
Output
|
|
<resample - upsample >
Ex)
Input
|
Fs=1000;
A=1.5;
B=1;
f1=50;
f2=100;
t=0:1/Fs:1;
ifreqs=4;
ti= 0:1/(Fs*ifreqs):1;
x=A*cos(2*pi*f1*t)+B*cos(2*pi*f2*t);
y=resample(x,ifreqs,1);
plot(ti(1:100),y(1:100),'ro','MarkerFaceColor',[1 0 0]);;hold on;
stem(t(1:25),x(1:25),'MarkerFaceColor',[0 0 1]);hold off; % plot original signal
legend('resampled data','original data');
|
Output
|
|
< resample - downsample >
Ex)
Input
|
Fs=1000;
A=1.5;
B=1;
f1=50;
f2=100;
t=0:1/Fs:1;
ifreqs=4;
ti= 0:1/(Fs*ifreqs):1;
x=A*cos(2*pi*f1*t)+B*cos(2*pi*f2*t);
y=resample(x,ifreqs,1);
y1=resample(y,1,ifreqs);
plot(ti(1:100),y(1:100),'ro','MarkerFaceColor',[1 0 0]);;hold on;
stem(t(1:25),y1(1:25),'MarkerFaceColor',[0 0 1]);hold off; % plot original signal
legend('original data','resampled data');
|
Output
|
|
< Filter >
Apply the Digital Filter described as below. As you see in the mathematical description, you can apply both FIR and IIR applying the tab values.
Ex)
Input
|
a = 1;
b = [0.2 0.2 0.2 0.2 0.2];
t = linspace(0,5*pi,100);
x = sin(t) .+ 0.5*rand(1,100);
y = filter(b,a,x);
plot(x,'r-',y,'b-');
|
Output
|
|
< fir1>
Case 1 : b = fir1(N,[pass_L,pass_R])
// Generate N tap Bandpass Filter with passband [pass_L,pass_R]. 0 < pass_L,pass_R < 1
Ex)
Input
|
b = fir1(48,[0.25 0.75]);
stem(b);
figure;
freqz(b,1,512);
|
Output
|
|
Case 2 : b = fir1(N,pass_R,'low')
// Generate N tap Lowpass Filter with passband [0,pass_R]. 0 < pass_R < 1
Ex)
Input
|
b = fir1(48,0.3,'low');
stem(b);
figure;
freqz(b,1,512);
|
Output
|
|
Case 3 : b = fir1(N,pass_L,'high')
// Generate N tap Highpass Filter with passband pass_L. 0 < pass_L < 1
Ex)
Input
|
b = fir1(48,0.7,'high');
stem(b);
figure;
freqz(b,1,512);
|
Output
|
|
Case 4 : b = fir1(N,[stop_L stop_R],'stop')
// Generate N tap Bandstop Filter with bandstop between stop_L and stop_R. 0 < stop_L, stop_R < 1
Ex)
Input
|
b = fir1(48,[0.45 0.55],'stop');
stem(b);
figure;
freqz(b,1,512);
|
Output
|
|
< fir2 >
Case 1 : b = fir2(N,[f1,f2,f3,...,fm],[mag1,mag2,mag3,...,magM])
// Generate N tap Filter with aribitrary passband definition specified by [f1,f2,f3,...,fm],[mag1,mag2,mag3,...,magM]
Ex)
Input
|
f=[0, 0.2, 0.2, 0.4, 0.4, 0.7, 0.7, 0.8, 0.8, 1];
m=[0, 0, 1, 1/2, 0, 0, 1, 1, 0, 0];
b=fir2(120,f,m);
figure;
plot(f,m,'r-');axis([0 1 -0.2 1.2]);
figure;
stem(b);axis([0,length(b),-0.5,0.5]);
figure;
freqz(b,1,512);
|
Output
|
|
< freqz >
Case 1 : freqz(b,a,NoOfFrequencyPoints)
// Produce frequency response plot (magnitude and phase vs frequency) with IIR filter coefficient vector b and a. (In case of FIR filter, set a to be 1).
Ex)
Input
|
b = fir1(48,[0.25 0.75]);
freqz(b,1,512);
|
Output
|
|
Case 2 : [h,w] = freqz(b,a,NoOfFrequencyPoints)
// return the numerical data of frequency response with IIR filter coefficient vector b and a. (In case of FIR filter, set a to be 1). h returns the magnitude values and w returns frequency point values.
Ex)
Input
|
b = fir1(48,[0.25 0.75]);
[h,w]=freqz(b,1,512);
plot(w,10 .* log(abs(h)),'r-');axis([0 pi,-100,0]);
|
Output
|
|
< RRC (Root Raised Cosine) Filter >
Ex) This example may look a little complicated but the most important parts are only two lines (rcosine(), rcosflt()) as marked in blue.
Input
|
sig = randi([0 1],1024,1); % This would generate random data of {0, 1}
sig = round(2*(sig - 0.5)); % This is to convert the data {0,1} into {-1,1}
upsample=8;
num = rcosine(1,upsample,'sqrt'); % Transfer function of filter
% This function performs two procedures in a single step. First it upsample the input data(sig) and filter it with rrc filter (num)
sig_rrc = rcosflt(sig,1,upsample,'filter',num);
% This is to perform fft. nextpow2() returns the number which is power of 2 and nearest to the number of input data.
% fftshift rotate the fft data so that the peak is located at the center of the spectrum.
Lsig = length(sig);
Lsig_rrc = length(sig_rrc);
NFFTsig = 2^nextpow2(Lsig);
NFFTsig_rrc = 2^nextpow2(Lsig_rrc);
sigfft=fftshift(fft(sig,NFFTsig)/Lsig);
sig_rrcfft=fftshift(fft(sig_rrc,NFFTsig_rrc)/Lsig_rrc);
% This is to convert the fft result into dB scale.
sigfft_dB = 20*log10(abs(sigfft));
sig_rrcfft_dB = 20*log10(abs(sig_rrcfft));
% this is for plotting the original data and filtered data and frequency spectrum
plotRange_tsig = 1:50;
plotRange_tsig_rrc = 1:(50*upsample);
subplot(2,2,1);plot(sig(plotRange_tsig));
axis([1 length(plotRange_tsig) -1.2 1.2]);title('original data');
subplot(2,2,2);plot(sigfft_dB);
axis([0 length(sigfft_dB) -100 0]);title('FFT for original data');
subplot(2,2,3);plot(sig_rrc(plotRange_tsig_rrc));
axis([1 length(plotRange_tsig_rrc) -1.2 1.2]);title('upsampled and filtered data');
subplot(2,2,4);plot(sig_rrcfft_dB);
axis([0 length(sig_rrcfft_dB) -100 0]);title('FFT for upsampled and filtered data');
|
Output
|
|
Ex) This example may look a little complicated but the most important parts are only two lines (rcosine(), rcosflt()) as marked in blue.
Input
|
data = randi([0 1],1024,1);
hMod = comm.QPSKModulator('BitInput',true)
hMod.PhaseOffset = 0;
sig = step(hMod, data);
upsample=8;
num = rcosine(1,upsample,'sqrt'); % Transfer function of filter
% This function performs two procedures in a single step. First it upsample the input data(sig) and filter it with rrc filter (num)
sig_rrc = rcosflt(sig,1,upsample,'filter',num);
% This is to perform fft. nextpow2() returns the number which is power of 2 and nearest to the number of input data.
% fftshift rotate the fft data so that the peak is located at the center of the spectrum.
Lsig = length(sig);
Lsig_rrc = length(sig_rrc);
NFFTsig = 2^nextpow2(Lsig);
NFFTsig_rrc = 2^nextpow2(Lsig_rrc);
sigfft=fftshift(fft(sig,NFFTsig)/Lsig);
sig_rrcfft=fftshift(fft(sig_rrc,NFFTsig_rrc)/Lsig_rrc);
% This is to convert the fft result into dB scale.
sigfft_dB = 20*log10(abs(sigfft));
sig_rrcfft_dB = 20*log10(abs(sig_rrcfft));
% this is for plotting the original data and filtered data and frequency spectrum
subplot(2,2,1);plot(real(sig),imag(sig));
axis([-1 1 -1 1]);title('original data');
subplot(2,2,2);plot(sigfft_dB);
axis([0 length(sigfft_dB) -100 0]);title('FFT for original data');
subplot(2,2,3);plot(real(sig_rrc),imag(sig_rrc));
axis([-1 1 -1 1]);title('upsampled and filtered data');
subplot(2,2,4);plot(sig_rrcfft_dB);
axis([0 length(sig_rrcfft_dB) -100 0]);title('FFT for upsampled and filtered data');
|
Output
|
|
Ex) This example may look a little complicated but the most important parts are only two lines (rcosine(), rcosflt()) as marked in blue.
Input
|
data = randi([0 1],1024,1);
hMod = comm.QPSKModulator('BitInput',true)
hMod.PhaseOffset = 0;
sig = step(hMod, data);
upsample=8;
num = rcosine(1,upsample,'sqrt'); % Transfer function of filter
% This function performs two procedures in a single step. First it upsample the input data(sig) and filter it with rrc filter (num)
sig_rrc = rcosflt(sig,1,upsample,'filter',num);
sig_rcv = rcosflt(sig_rrc,1,upsample,'Fs/filter',num);
% This is to perform fft. nextpow2() returns the number which is power of 2 and nearest to the number of input data.
% fftshift rotate the fft data so that the peak is located at the center of the spectrum.
Lsig = length(sig);
Lsig_rrc = length(sig_rrc);
Lsig_rcv = length(sig_rcv);
NFFTsig = 2^nextpow2(Lsig);
NFFTsig_rrc = 2^nextpow2(Lsig_rrc);
NFFTsig_rcv = 2^nextpow2(Lsig_rcv);
sigfft=fftshift(fft(sig,NFFTsig)/Lsig);
sig_rrcfft=fftshift(fft(sig_rrc,NFFTsig_rrc)/Lsig_rrc);
sig_rcvfft=fftshift(fft(sig_rcv,NFFTsig_rcv)/Lsig_rcv);
% This is to convert the fft result into dB scale.
sigfft_dB = 20*log10(abs(sigfft));
sig_rrcfft_dB = 20*log10(abs(sig_rrcfft));
sig_rcvfft_dB = 20*log10(abs(sig_rcvfft));
% this is for plotting the original data and filtered data and frequency spectrum
subplot(3,2,1);plot(real(sig),imag(sig));
axis([-1 1 -1 1]);title('original data');
subplot(3,2,2);plot(sigfft_dB);
axis([0 length(sigfft_dB) -100 0]);title('FFT for original data');
subplot(3,2,3);plot(real(sig_rrc),imag(sig_rrc));
axis([-1 1 -1 1]);title('upsampled and filtered data');
subplot(3,2,4);plot(sig_rrcfft_dB);
axis([0 length(sig_rrcfft_dB) -100 0]);title('FFT for upsampled and filtered data');
subplot(3,2,5);plot(real(sig_rcv),imag(sig_rcv));
axis([-2 2 -2 2]);title('recieved data');
subplot(3,2,6);plot(sig_rcvfft_dB);
axis([0 length(sig_rcvfft_dB) -100 0]);title('FFT for recieved data');
|
Output
|
|
|
|