Sunday, January 31, 2010

DSP: Using a basic digital filter: Stock market example


Digital filters are very powerful tools in DSP. In this example we will use a digital filter to help with stock market analysis. In the figure above, you can see the raw data plotted with a 5 day moving average and a 20 day moving average. These two moving average lines can be used for buy/sell predictions.

To calculate the 20 day moving average, the window size and filter coefficients must be established. The following Matlab code sets the window size and calculates the filter results. The five day filter is implemented using the same code, with a windowSize = 5. The matrix x in the example contains the raw stock data. For more info on the filt function: http://www.mathworks.com/access/helpdesk/help/techdoc/ref/filter.html

windowSize = 20; % Twenty day average
% The following filter applies the 1/windowSize coefficients
twentyDay = filter(ones(1,windowSize)/windowSize,1,x);


Thursday, January 28, 2010

Use of Windowing Techniques on 3G CDMA Signal


Here is an example of how windowing (blackman in this example) can help reduce sidebands. This technique is useful in making spectral measurements such as Adjacent Channel Leakage Ratio. By minimizing the sidebands generated by an arbitrary waveform generator, for example, the true device performance can be more accurately measured. For more information on windowing, see the following article from National Instruments: http://zone.ni.com/devzone/cda/tut/p/id/4844

Tuesday, January 26, 2010

Coherent Noise Averaging




% James Eastham
% Member, IEEE
% Revision: R1
% Contact: james.eastham@ieee.org
%
% Description: This script shows the affect of coherent averaging to
% reduce random noise.

%% Affect of Additive Noise
close all; % close all plots
clear all; % clear all variables
fs = 500; % sample freq
Ts=1/fs; % sample period
n=0:1*fs; % time window (fs=samples per second), 1 second
tn=n*Ts; % sample index
x=cos(2*pi*50*tn); % 50Hz sinusoid
k=.5; %noise amplitude
n=randn(size(x)); % random noise
s=x+(k*n); %signal + random noise with amplitude k
% Time domain plot of x(t)=cos(2*pi*50*tn)
figure('Color',[1 1 1]);
plot(tn,s,'g','LineWidth',2);box off;hold on;
title('Sinusoid w/ Additive White Noise');
xlabel('Time(s)');
ylabel('Amplitude');
% Frequency Domain Plot
figure('Color',[1 1 1]);
[X,f]=ComputeSpectrum(s,fs,2^12);
orig_X=X;
h=plot(f,20*log(X));
box off;
grid on;
set(h,'Linewidth',1);
xlabel('Frequency(Hz)');
ylabel('Amplitude(dB)');
title(['Frequency Spectrum w/ Additive Noise k=.5']);
axis([0 250 -50 120]);

%% Example of Noise Averaging, 3 Coherent Averages
fs = 500; % sample freq
Ts=1/fs; % sample period
n=0:1*fs; % time window (fs=samples per second), 1 second
tn=n*Ts; % sample index
x=cos(2*pi*50*tn);
k=.2; %Noise amplitude
n1=randn(size(x));
n2=randn(size(x));
n3=randn(size(x));
s1=x+(k*n1);
s2=x+(k*n2);
s3=x+(k*n3);
average_s=((s1+s2+s3)/3);
% Time domain plot of x(t)=cos(2*pi*50*tn)
figure('Color',[1 1 1]);
h=plot(tn,average_s,'g','LineWidth',2);box off;hold on;
title('Sinusoid w/ Additive White Noise: 3 Averages');
xlabel('Time(s)');
ylabel('Amplitude');
% Frequency Domain Plot
figure('Color',[1 1 1]);
[X,f]=ComputeSpectrum(average_s,fs,2^12);
plot(f,20*log(X),f,20*log(orig_X),'--r');
box off;
grid on;
set(h,'Linewidth',1);
xlabel('Frequency(Hz)');
ylabel('Amplitude(dB)');
legend('3 Averages','Original');
title(['Frequency Spectrum w/ Noise Averaging: 3 Averages']);
axis([0 250 -50 120]);

%% Example of Noise Averaging, 10 Coherent Averages
fs = 500; % sample freq
Ts=1/fs; % sample period
n=0:1*fs; % time window (fs=samples per second), 1 second
tn=n*Ts; % sample index
x=cos(2*pi*50*tn);
k=.2; %Noise amplitude
n1=randn(size(x));
n2=randn(size(x));
n3=randn(size(x));
n4=randn(size(x));
n5=randn(size(x));
n6=randn(size(x));
n7=randn(size(x));
n8=randn(size(x));
n9=randn(size(x));
n10=randn(size(x));
s1=x+(k*n1);
s2=x+(k*n2);
s3=x+(k*n3);
s4=x+(k*n1);
s5=x+(k*n2);
s6=x+(k*n3);
s7=x+(k*n1);
s8=x+(k*n2);
s9=x+(k*n3);
s10=x+(k*n1);
average_s=((s1+s2+s3+s4+s5+s6+s7+s8+s9+s10)/10);
% Time domain plot of x(t)=cos(2*pi*50*tn)
figure('Color',[1 1 1]);
h=plot(tn,average_s,'g','LineWidth',2);box off;hold on;plot(tn,x,'--r');
title('Sinusoid w/ Additive White Noise: 10 Averages');
xlabel('Time(s)');
ylabel('Amplitude');
legend('composite signal','averaged');
% Frequency Domain Plot
figure('Color',[1 1 1]);
[X,f]=ComputeSpectrum(average_s,fs,2^12);
plot(f,20*log(X),f,20*log(orig_X),'--r');
box off;
grid on;
set(h,'Color',[0 0.8 1]);
set(h,'Linewidth',1);
xlabel('Frequency(Hz)');
ylabel('Amplitude(dB)');
title(['Frequency Spectrum w/ Noise Averaging: 10 Averages']);
legend('10 Averages','Original');
axis([0 250 -50 120]);

Friday, January 22, 2010

DSP Tutorial #5: Frequency Analysis


This example computes the spectrum of a complex signal using the ComputeSpectrum function.

Couple important things to notice:
- Our range of frequencies on the x-axis is Fs/2
- As Fs increases, the resolution in the frequency domain is much better
- As our FFT bins (2^N) increases, the more bins we have (we can still resolve the freq content provided we meet the sample theorem, we have more bins - basically a finer plot)
- As our windows length increases, the better resolution we have


close all;
clear all;
fs=500;
Ts=1/fs;
n=0:1*fs;
tn=n*Ts;
x=5+cos(2*pi*200*tn)+10+cos(2*pi*100*tn);
% Time Domain Plot
figure('Color',[1 1 1]);
plot(tn,x);box off;
title(['Time Domain Signal']);
[X,f]=ComputeSpectrum(x,fs,2^14);
figure('Color',[1 1 1]);
plot(f,20*log(X));
title(['Freq Domain']);
axis('tight');

The Compute Spectrum Function

This function uses the FFT algorithm to compute the spectrum of an input with a number of points:

function [X,f]=ComputeSpectrum(x,fs,N);
% ComputeSpectrum Compute the Spectrum
%
% x input signal
% fs Sampling frequency of the input signal
% N number of points to evaluate spectrum (integer number of 2)
%
% X Absolute value of the spectrum (capital X for spectrum)
% f Frequency axis of the spectrum
%
% This function uses the FFT algorithm to compute the spectrum of an
% input signal with a number of points (N) N-point FFT where N must be
% a power of 2 (e.g. 2^10)
%
%Example: Compute the spectrum of sinusoid at 3Hz, fs=36Hz
%
% x=cos(2*pi*3*[0:12*3]*1/(12*3)); 12 times max freq
% [X,f]=ComputeSpectrum(x,12*3,2^12);
%
%Revision:0.0.1 jeastham
%
X=abs(fft(x,N));
f=(fs/N)*[0:N-1]; %make the axis freq
X=X(1:end/2);
f=f(1:end/2);

Monday, January 18, 2010

DSP Tutorial #4: Plotting Complex Signals

In this example we will plot a more complex time domain signal in Matlab:



clear all; %clear all variables
fs = 50; % sample freq
Ts=1/fs; % sample period
n=0:1*fs; % time window (fs=samples per second), 1 second
tn=n*Ts; % sample index
sig_a=cos(2*pi*5*tn);
sig_b=1.5+cos(2*pi*20*tn);
x=sig_a+sig_b;
% Time domain plot of x(t)=cos(2*pi*5*tn)
figure('Color',[1 1 1]);
h=plot(tn,x,'g','LineWidth',3);box off;hold on;
plot(tn,sig_a,'--r','LineWidth',2);
plot(tn,sig_b,'--k','LineWidth',2);
title('Example of Composite Signal: 2 Signals Added');
legend('Composite','Signal 1','Signal 2');
xlabel('Time(s)');
ylabel('Amplitude');
figure('Color',[1 1 1]);
plot(tn,sig_a);
title('Stem Plot Example: Fs = 50 (50 Points)');
axis('tight');

Sum of Sinusoids Triangle Approximation


Triangle Approximation: Sum of Sinusoids Example

The following Matlab code provides an example of the fact that any signal can constructed using a sum of sinusoids.



clear all; %clear all variables
fs = 50; % sample freq
Ts=1/fs; % sample period
n=0:1*fs; % time window (fs=samples per second), 1 second
tn=n*Ts; % sample index
wave_1=sin(2*pi*5*tn);
wave_2=sin(-1/9*sin(6*pi*5*tn));
wave_3=sin(1/25*sin(10*pi*5*tn));
wave_4=sin(-1/49*sin(14*pi*5*tn));
wave_5=sin(1/81*sin(18*pi*5*tn));
wave_6=sin(-1/121*sin(22*pi*5*tn));
wave_7=sin(1/169*sin(26*pi*5*tn));
x = (8/(pi^2))*(wave_1+wave_2+wave_3+wave_4+wave_5+wave_6+wave_7);
% Time domain plot of x(t)=cos(2*pi*5*tn)
figure('Color',[1 1 1]);
h=plot(tn,wave_1,'--r',tn,wave_2,'--k',tn,x,'g','LineWidth',2);box off;hold on;
xlabel('Time(s)');
ylabel('Amplitude');
title('Square Wave Approximation');
axis([0 1 -1.5 1.5]);