Wednesday, May 7, 2014

Implementing filters in Matlab using ellipord and filtfilt


We can create a custom function in Matlab using the built in ellipord and filtfilt functions included in the signal processing toolbox.

Here is a lowpass filter with explanation in the comments. Once we create this lowpass function, we can call this function and simply pass the function our input signal, sampling frequency, and cut-off frequency.

function [y,N,B,A]= lowpass(x,fs,fc);
% Description: This function uses an elliptical filter
% to filter out high frequencies.
%
% [y,N,B,A] = lowpass(x,fs,fc)
%
% Variable Definition:
% x: Input signal
% fs: Input signal sampling frequency in Hz
% fc: cutoff frequency
%
% y: Filtered signal (output)
% N: Filter order
% B: Feedforward (FIR) coefficients
% A: Feedback (IIR) coefficients
Rp = 0.5;   % Ripple in dB (passband)
Rs = 20;    % Attenuation
mf = fs/2;  % Maximum frequency
Wp = fc/mf; % Cuttoff frequency
Ws = 1.2 * Wp; %Stopband frequency
% First, determine the minimum filter order to meet specifics
[N,Wn] = ellipord(Ws,Wp,Rp,Rs);
% Next, Design the filter (find the coefficients)
[B,A] = ellip(N,Rp,Rs,Wn,'low');
% Now, apply the filter (filter the signal)
y = filtfilt(B,A,x);
% Plot the filter
figure('Color',[1 1 1]);
freqz(B,A,2^14,fs);

Monday, November 11, 2013

Automatic Sorting in Excel without a Macro (Sort Chart Data)

I really wanted to have a chart on my sheet which was sorted from high to low. Imagine if you wanted to create a bar chart which showed net income by company and wanted the chart to be sorted high to low. One solution is to just manually sort the data. Another option is to create a macro which copies the data somewhere (values only) and sorts the data, then chart the sorted data. These options all require the user to 'do something' to get the sorted chart.

I came up with a new approach using just the MAX feature, VLOOKUP, and a simple IF. The original data is in the left columns C3:C8. I just take the MAX of the group, then use a simple IF statement to assign a "One", "Two", etc, then just a VLOOKUP to get the order back at the bottom of the sheet. The IF statement is just =IF(C4=$C$9,"One",C4).

The reason I spell out the rank order is because if you don't do this and have negative numbers the logic will get screwed up as rank "5" will be greater than 99 for example and mess up your order.

Let me know if you have questions, the chart is then done on the table below. Pretty simple.

Monday, September 9, 2013

Autocreate Minitab Scripts in C++

The following code will read in a CSV file containing row and reference line data. The program will then create a Minitab script as the output file. This file can then be used directly in Minitab to auto-generate hundreds of plots. This approach can work for auto-generating / automating any scripts in Minitab.

// main_histo-limits.cpp
// James Eastham
// Version 1.0
// Purpose: The purpose of this program is to read in a set of test names and
// limits. Produce minitab script code which can be used to generate a large
// number of histograms with correct reference limit lines.
// Input: The input to this program will be a CSV file containing all the test
// names and LSL/USL for which histograms will be created.
// The input file is hard coded as limits.csv
// Output: The output will be a text file (histolimits.txt) which can be pasted
// directly into minitab and ran as a script
#include
#include
#include ;
#include ;The

using namespace std;

// int argc, char *argv[] are only needed for processing command line
// arguments. This program does not use command line arguments, so
// main() would be fine here.

int main()
{
    string testname$;           // Variable for test name read from CSV
    string lsl$;                // Varialbe for lower spec limit from CSV
    string usl$;                // Variable for upper spec limit from CSV
    // Input file name is hardcoded below as testnames.txt
    ifstream testnames ("limits.csv");
    ofstream myfile;
    // Output script file is hardcoded below as histoscript.txt
    myfile.open ("histolimits.txt");
    if (testnames.is_open())
    {
       while (testnames.good() )
             {
             getline (testnames,testname$,',');  // Read test name from file
             getline (testnames,lsl$,',');       // Read lsl from file
             getline (testnames,usl$);           // Read usl from file
             // The following cout lines simply print the info read from CSV
             // to the screen for debug purposes
             cout << testname$;  
             cout << ',';
             cout << lsl$;
             cout << ',';
             cout << usl$;
             cout << endl;
             if (testname$.length()>30)
                {
                cout << "***WARNING LONG TEST NAME, Exceeds 30 Charcaters***\n";
                }
    // The following 'myfile' lines output text to the output file
    // 'histolimits.txt'. This txt file is then used in Minitab as
    // a script. The testname$, lsl$, usl$ are used in this Minitab script
    // for naming each histogram/plot and reference/limit lines.
    // This text file can be pasted directly into Minitab.
   
    myfile << "Histogram '";
    myfile << testname$;
    myfile << "';\n";
    myfile << "Scale 1;\n";
    myfile << "TFont \"Arial Narrow\";\n";
    myfile << "PSize 8;\n";
    myfile << "Scale 2;\n";
    myfile << "TFont \"Arial Narrow\";\n";
    myfile << "PSize 8;\n";
    myfile << "AxLabel 1;\n";
    myfile << "TFont \"Arial Narrow\";\n";
    myfile << "AxLabel 2;\n";
    myfile << "TFont \"Arial Narrow\";\n";
    myfile << "PSize 8;\n";
    myfile << "NoBold;\n";
    myfile << "Table;\n";
    myfile << "TFont \"Arial Narrow\";\n";
    myfile << "Section 1;\n";
    myfile << "NoLegend;\n";
    myfile << "Density;\n";
    myfile << "Bar 'part_id';\n";
    myfile << "Color 3;\n";
    myfile << "Distribution 'part_id';\n";
    myfile << "Normal;\n";
    myfile << "Grid 1;\n";
    myfile << "Reference 1 ";
    myfile << lsl$;
    myfile << ";\n";
    myfile << "Type 3;\n";
    myfile << "Color 2;\n";
    myfile << "Size 1;\n";
    myfile << "MODEL 1;\n";
    myfile << "TFont \"Arial Narrow\";\n";
    myfile << "PSize 8;\n";
    myfile << "Reference 1 ";
    myfile << usl$;
    myfile << ";\n";
    myfile << "Type 3;\n";
    myfile << "Color 114;\n";
    myfile << "Size 1;\n";
    myfile << "MODEL 1;\n";
    myfile << "Title \"TM6PP6924 : ";
    // The following two remarks are not needed
    //myfile << "Title &\n";
    //myfile << "\"";
    myfile << testname$;
    myfile << "\";\n";
    myfile << "TFont \"Arial Narrow\";\n";
    myfile << "PSize 12;\n";
    myfile << "NoBold;\n";
    myfile << "SubTitle;\n";
    myfile << "TFont \"Arial Narrow\";\n";
    myfile << "PSize 8;\n";
    myfile << "StDist;\n";
    myfile << "Footnote;\n";
    myfile << "FPanel;\n";
    myfile << "NoDTitle;\n";
    myfile << "NoDSubtitle.\n";
     
}
testnames.close();   // close the input file
}
else cout << "Unable to open file";
   
    myfile.close();  // close the output file
    system("PAUSE");
    return EXIT_SUCCESS;
}

Monday, October 15, 2012

Matlab Interpolation

Here is a quick example of linear interpolation.




%% Linear Interpolation Example
% James Eastham, Member IEEE
% 10/15/2012
close all;
clear all;
x = 0:100; %Original Data Points
y = sin(x); % sin of original data points
xi = 0:.5:100; % points to interpolate
yi = interp1(x,y,xi); % interpolated vector
figure('Color',[1 1 1]);
subplot(3,1,1);
plot(x,y,'o-g',xi,yi,'.-b')
title('Interpolation Example 1 Points');
xlabel('point');
ylabel('sin(x)');
x = 0:100; %Original Data Points
y = sin(x); % sin of original data points
xi = 0:.25:100; % points to interpolate
yi = interp1(x,y,xi); % interpolated vector
subplot(3,1,2);
plot(x,y,'o',xi,yi,'.-r')
title('Interpolation Example 3 Points');
xlabel('point');
ylabel('sin(x)');
x = 0:100; %Original Data Points
y = sin(x); % sin of original data points
xi = 0:.1:100; % points to interpolate
yi = interp1(x,y,xi); % interpolated vector
subplot(3,1,3);
plot(x,y,'o',xi,yi,'.-g')
title('Interpolation Example 9 Points');
xlabel('point');
ylabel('sin(x)');

Tuesday, May 1, 2012

Terrain Profile & Path Loss



The site "Hey, what's that?" (http://www.heywhatsthat.com/) provides a great free tool for terrain profiles, coverage, and basic path calculations.

The Path Profiler is located at:
http://www.heywhatsthat.com/profiler-0904.html


Friday, April 20, 2012

Fourier Additivity Property



   



 The plot above shows two signals being added in the time domain and the spectrum of each time domain signal. The plot below shows the spectrum of each signal added in the frequency domain.




% Here we demonstrate the concept of additivity, that is, adding two or
% more signals in one domain results in the corresponding signals being
% added in the other domain.
close all;
    fs=125;   %sample frequency = 125Hz
    Ts=1/fs;  %sample period
    n=0:100*fs; %sample index (fs=samples per second), w=window lenght
    tn=n*Ts;
    x=cos(2*pi*5*tn);
% Time domain plot of x(t)=cos(2*pi*5*tn)
    figure('Color',[1 1 1]);
    subplot(3,2,1);
    h=plot(tn,x);box off;grid on;
    xlabel('Time(s)');
    ylabel('Amplitude');
    title(['Time Domain: x1(t)=cos(2*pi*5*tn)']);
    axis([0 1 -1.2 1.2]);
% Frequency domain plot of x(t)=cos(2*pi*5*tn)
    subplot(3,2,2);
    [X,f]=ComputeSpectrum(x,fs,2^14);
    sig1=X;
    h=plot(f,20*log(X));
    box off; grid on;
    set(h,'Linewidth',1);
    xlabel('Frequency(Hz)');
    ylabel('Amplitude(dB)');
    title(['Frequency Domain: x1(t)=cos(2*pi*5*tn)']);
    axis([4.4 5.8 -100 200]);
   
    x=cos(2*pi*10*tn);
  % Time domain plot of x(t)=cos(2*pi*10*tn)
    subplot(3,2,3);
    h=plot(tn,x);box off;grid on;
    xlabel('Time(s)');
    ylabel('Amplitude');
    title(['Time Domain: x2(t)=cos(2*pi*10*tn)']);
    axis([0 1 -1.2 1.2]);
% Frequency domain plot of x(t)=cos(2*pi*10*tn)
    subplot(3,2,4);
    [X,f]=ComputeSpectrum(x,fs,2^14);
    h=plot(f,20*log(X));
    sig2=X;
    box off; grid on;
    set(h,'Linewidth',1);
    xlabel('Frequency(Hz)');
    ylabel('Amplitude(dB)');
    title(['Frequency Domain: x2(t)=cos(2*pi*10*tn)']);
    axis([5 15 -100 250]);

    x=(cos(2*pi*5*tn)) + (cos(2*pi*10*tn)); % add two signal in time domain
% Time domain plot of x(t)=cos(2*pi*5*tn)+cos(2*pi*10tn)
    subplot(3,2,5);
    h=plot(tn,x);box off;grid on;
    xlabel('Time(s)');
    ylabel('Amplitude');
    title(['Time Domain: x1(t) + x2(t)']);
    axis([0 1 -3 3]);
% Frequency domain plot of x(t)=cos(2*pi*5*tn)+cos(2*pi*10tn)
    subplot(3,2,6);
    [X,f]=ComputeSpectrum(x,fs,2^14);
    h=plot(f,20*log(X));
    box off; grid on;
    set(h,'Linewidth',1);
    xlabel('Frequency(Hz)');
    ylabel('Amplitude(dB)');
    title(['Frequency Domain: x1(t) + x2(t)']);
    axis([0 30 -100 250]);
 % Now, let's simply add the specturm of each signa.
    figure('Color',[1 1 1]);
    subplot(3,1,1);
    h=plot(f,20*log(sig1));
    box off; grid on;
    set(h,'Linewidth',1);
    xlabel('Frequency(Hz)');
    ylabel('Amplitude(dB)');
    title(['Frequency Domain: x1(t) + x2(t)']);
    axis([0 30 -100 250]);
   
    subplot(3,1,2);
    h=plot(f,20*log(sig2));
    box off; grid on;
    set(h,'Linewidth',1);
    xlabel('Frequency(Hz)');
    ylabel('Amplitude(dB)');
    title(['Frequency Domain: x1(t) + x2(t)']);
    axis([0 30 -100 250]);

    subplot(3,1,3);
    h=plot(f,20*log(sig1+sig2));
    box off; grid on;
    set(h,'Linewidth',1);
    xlabel('Frequency(Hz)');
    ylabel('Amplitude(dB)');
    title(['Frequency Domain: x1(t) + x2(t)']);
    axis([0 30 -100 250]);

Sunday, February 5, 2012

Normal Distributions in MatLab

%% Normal Distributions in Matlab
% James Eastham
% Member, IEEE
% Created on: 2/5/2012
%
% The following code creates multiple normal distributions

figure('Color',[1 1 1]);
x = 17:0.1:28; %Plot range

pop_average = 21.98; %population mean
sigma = .4; %population sigma
normal_dist = normpdf(x, pop_average, sigma);
plot(x, normal_dist,'LineWidth',3);hold on;

pop_average = 22.89;
sigma = .244;
normal_dist = normpdf(x, pop_average, sigma);
plot(x, normal_dist,'--r','LineWidth',3);hold on;

pop_average = 21.84;
sigma = .78;
normal_dist = normpdf(x, pop_average, sigma);
plot(x, normal_dist,'-.g','LineWidth',3);hold on;

pop_average = 22.24;
sigma = .91;
normal_dist = normpdf(x, pop_average, sigma);
plot(x, normal_dist,'-ko','LineWidth',3);hold on;

legend('build 1','build 2','build 3','overall')
grid on;
title('Design Build Performance');
xlabel('Gain(dB)');