freq_plot.m

%a function which plots the magnitude spectrum of a signal from DC (0Hz) up to
%max_freq.
%
%   example: 
%           fs = 10000;
%           t = 0:1/fs: 4; 
%           sig = 3*cos(2*pi*100*t) + 3*cos(2*pi*20*t); 
%           freq_plot(sig, 120, fs)
%
% by david.dorran@gmail.com, Feb 2012
function freq_plot(sig, max_freq, fs)
    if(max_freq>fs/2)
        error('Max freq must be less than Nyquist frequency');
    end
    X = abs(fft(sig));
    N = length(sig);
    
    f = 0:fs/length(sig):fs;
    num_bins = find(f> max_freq);
    X = X(1:num_bins(1));
    f = f(1:num_bins(1));
    plot(f, X/(N)*2+0.0001)
    xlabel('frequency (Hz)')
    ylabel('Magnitude');

Advertisement

Basic Matlab GUI

Here is the code I used to create a basic matlab GUI without using GUIDE in the following youtube videos:
First Video: GUI Layout http://youtu.be/kDcn83becS0
Second Video: Object handles and Object Manipulation http://youtu.be/XGRkqgDwWRA
Third Video: Callback Functions http://youtu.be/JT8gAbY2Ooc
Fourth Video: Finishing Touches http://youtu.be/VCyKQAMUpmU

% Basic Graphical User Interface (GUI) without using GUIDE
% Available at https://dadorran.wordpress.com

% There are three basic areas to understand:
%   1. Layout    (how to position objects on the GUI)
%   2. Handles to Objects (used to modify object properties)
%   3. Callback functions (used by User Interface Objects)


%create a figure to house the GUI
figure

                
%create an annotation object 
ellipse_position = [0.4 0.6 0.1 0.2];
ellipse_h = annotation('ellipse',ellipse_position,...
                    'facecolor', [1 0 0]);
                
%create an editable textbox object
edit_box_h = uicontrol('style','edit',...
                    'units', 'normalized',...
                    'position', [0.3 0.4 0.4 0.1]);

%create a "push button" user interface (UI) object
but_h = uicontrol('style', 'pushbutton',...
                    'string', 'Update Color',...
                    'units', 'normalized',...
                    'position', [0.3 0 0.4 0.2],...
                    'callback', {@eg_fun,edit_box_h, ellipse_h });
            
%Slider object to control ellipse size
uicontrol('style','Slider',...        
            'Min',0.5,'Max',2,'Value',1,...
            'units','normalized',...
            'position',[0.1    0.2    0.08    0.25],...
            'callback',{@change_size,ellipse_h,ellipse_position });
        
uicontrol('Style','text',...
            'units','normalized',...
            'position',[0    0.45    0.2    0.1],...
            'String','Ellipse Size')

%eg_fun code ---------------------------------------------------
%copy paste this code into a file called eg_fun.m
function eg_fun(object_handle, event)
    disp('hi')

%updated eg_fun used to demonstrate passing variables
%copy paste this code into a file called eg_fun.m
function eg_fun(object_handle, event, edit_handle, ellipse_handle)
    str_entered = get(edit_handle, 'string');
    
    if strcmp(str_entered, 'red')
        col_val = [1 0 0];
    elseif strcmp(str_entered, 'green')
        col_val = [0 1 0];
    elseif strcmp(str_entered, 'blue')
         col_val = [0 0 1];
    else
        col_val = [0 0  0];
    end
    set(ellipse_handle, 'facecolor', col_val);
                
%change_size code --------------------------------------------------
%copy paste this code into a file called change_size.m
function  change_size(objHandel, evt, annotation_handle, orig_pos)
    slider_value = get(objHandel,'Value');
    new_pos = orig_pos;
    new_pos(3:4) = orig_pos(3:4)*slider_value;
    set(annotation_handle, 'position', new_pos)

DFT Windowing

This is code I used in a youtube video on DFT windowing – see http://youtu.be/PYd-pzaZpYE

% Tutorial on why DFT windowing works.
% See http://youtu.be/V_dxWuWw8yM for intro video
% Code available at https://dadorran.wordpress.com. 

fs = 1000;
T = 1/fs;
N = 1000; % Create a signal with N samples = 1 second
n = 0:N-1; % sample numbers
t = n*T; % sampling times over a one second duration.
x = cos(2*pi*6.5*t);

x_zpad = [x zeros(1, 99000)]; %zero pad
X_ZPAD_mags= abs(fft(x_zpad));

subplot(2,1,1)
plot(n,x);
xlabel('Samples');
ylabel('Amplitude')
title('Time-Domain Signal to be analysed');

subplot(2,1,2)
bins = 0:length(x_zpad)-1;
plot(bins(1:1500), X_ZPAD_mags(1:1500))
xlabel('Bins');
ylabel('Magnitude')
title('First 1500 Bins of Magnitude Spectrum of Zero-Padded Signal');

% Plot magnitude spectrum against frequency in Hertz
fax = [0:100000-1]*fs/100000;
plot(fax(1:1500), X_ZPAD_mags(1:1500));
title('Magnitude Spectrum of Zero-Padded Signal');
xlabel('Frequency (Hertz)');
ylabel('Magnitude')

%window using a hann window function
han_win = hanning(N)';
x_zpad_win = [x.*han_win zeros(1, 99000)]; %zero pad
X_ZPAD_WIN_mags= abs(fft(x_zpad_win));
hold on
plot(fax(1:1500), X_ZPAD_WIN_mags(1:1500),'r');


%The zero padded signal can be considered as being a sinusoid
% which is multiplied by a rectangular winddow. The folloiwing
% code illustrates this.
Ndisp = 5000; %number of time-domain samples to display
subplot(3,1,1)
plot(x_zpad(1:Ndisp))
title(['Zero Padded Signal (first ' num2str(Ndisp) ' samples)'])
ylabel('Amplitude')
xlabel('Samples');
t = [0:100000-1]*T; % sampling times for 100000 samples.
sinusoid = cos(2*pi*6.5*t);
subplot(3,1,2)
plot(sinusoid(1:Ndisp))
title(['Sinusoidal Signal (first ' num2str(Ndisp) ' samples)'])
ylabel('Amplitude')
xlabel('Samples');
subplot(3,1,3)
win = [ones(1, N) zeros(1, 99000)];
%win = [hanning(N)' zeros(1, 99000)];
plot(win(1:Ndisp))
title(['Rectangular Window (first ' num2str(Ndisp) ' samples)'])
ylabel('Amplitude')
xlabel('Samples');
ylim([-0.5 1.5])

%When you multiply signals in the time-domain it is equivalent to 
% convolution in the frequency domain. Let's take a look at each of 
% these three waveforms in the frequency domain
%update fax to  show negative frequencies
fax = [0:100000-1]*fs/100000;
fax(end-49998:end) = fax(end-49998:end)-1000;
fax(50001) = NaN;
subplot(3,1,1)
plot(fax,abs(fft(x_zpad)))
title(['Double Sided Magnitude Spectrum of  Zero Padded Signal'])
ylabel('Magnitude')
xlabel('Frequency (Hz)');

subplot(3,1,2)
plot(fax,abs(fft(sinusoid)))
title(['Double Sided Magnitude Spectrum of  Sinusoidal Signal'])
ylabel('Magnitude')
xlabel('Frequency (Hz)');

subplot(3,1,3)
plot(fax, abs(fft(win)))
title(['Double Sided Magnitude Spectrum of' ...
    ' Rectangular Window function '])
ylabel('Magnitude')
xlabel('Frequency (Hz)');

% Determine the magnitude spectrums with less zero padding
% Plot magnitude spectrum against frequency in Hz rather than bins
% Observe that the lower res spectrums are a sampled version of the
% higher res ones
X_mags4000 = abs(fft([x zeros(1, 3000)]));
X_mags2000 = abs(fft([x zeros(1, 1000)]));
X_mags1000 = abs(fft(x));

figure
fax = [0:4000-1]*fs/4000;
plot(fax(1:60), X_mags4000(1:60));
hold on
plot(fax(1:60), X_mags4000(1:60),'.');
set(gca,'XTick', fax(1:4:60))
xlabel('Frequency in Hertz')
ylabel('Magnitude')

fax = [0:2000-1]*fs/2000;
plot(fax(1:30), X_mags2000(1:30),'rx');

fax = [0:1000-1]*fs/1000;
plot(fax(1:15), X_mags1000(1:15),'go');

fax = [0:100000-1]*fs/100000;
plot(fax(1:1500), X_ZPAD_mags(1:1500),'k');

%Each of the other magnitude spectrums are a subset or sampled
% version of the high resolution one.
%It should be appreciated from these examples that all DFT's can
% be thought of as being a sampled version of the DTFT (which 
% produces a continuous spectrum that could be sampled at any 
% frequency)
%Understanding the effects of windowing therefore requires an
% undertstanding of what happening for the high res case.



% Comparison of some common window functions
figure
N = 1000;
zpad_len = 99000;
L = N+zpad_len;
zpad = zeros(1,zpad_len);
hamming_freq = fftshift(abs(fft([hamming(N)', zpad])))/(N);
plot([-5 :.01: 5],hamming_freq(L/2-500:L/2+500),'r')
hold on
rect_freq = fftshift(abs(fft([rectwin(N)', zpad])))/(N);
plot([-5 :.01: 5],rect_freq(L/2-500:L/2+500),'g')
hann_freq = fftshift(abs(fft([hanning(N)', zpad])))/(N);
plot([-5 :.01: 5],hann_freq(L/2-500:L/2+500),'c')
blackman_freq = fftshift(abs(fft([blackman(N)', zpad])))/(N);
plot([-5 :.01: 5],blackman_freq(L/2-500:L/2+500),'k')
legend('Hamming', 'Rectangular','Hann','Blackman')
ylabel('Magnitude');
xlabel('Bin width (non zero padded)');