DFT Basis Function Plots

Just posting this as I had request for code I used in creating a youtube video on the DFT

%This a script I used when creating the images for the DFT explained youtube video http://youtu.be/B2iUDBZzBpY
%
% It shows the the result of multiplying a signal being analysed by a set of sinusoidal basis functions (which is how the DFT works)
% The figure produced was edit further before using it in the video
close all;
N = 8; %signals are only 8 samples long
n = 0:N-1;

%specify the magnitudes and phases of the sinusoids (integer cyles over N samples)
mags = [ 0 0 5 0 0 0 ]; %change this to include more sinusoids i.e. change 0 values to non zero
phases = [0 0 0 0 0 0]; %phase values for each of the sinusoids

y_range = [-5 5];
interp_val = 100;
n_vis = 0 : (N)*interp_val-1;
x_vis = zeros(1, (N)*interp_val);
x = zeros(1, N);

for k = 0 : N/2
    x = x + mags(k+1)*cos(2*pi*(k)*n/N + phases(k+1)); %samples of signal
    x_vis = x_vis + mags(k+1)*cos(2*pi*(k)*n_vis/(N*interp_val) + phases(k+1)); %interpolated sinusoid to aid with visualisation
end

num_rows = N/2+2 +N/2; %number of rows in the subplot
subplot(num_rows,2,1)
plot(x_vis);
hold on
plot(n*interp_val, x,'r.')
hold off

set(gca, 'XTick', 0:interp_val:(N)*interp_val)
set(gca, 'XTickLabel', '')
set(gca, 'Box', 'off')
    set(gca, 'YTick', y_range)
        ylim(y_range)
%for each of the basis functions
for k = 0 : N/2
    basis_function = cos(2*pi*k*n/N); %samples of the cosine basis function
    basis_function_vis = cos(2*pi*k*n_vis/(N*interp_val)); %an interpolated version to aid with visualisation
    
    mult_sig = basis_function.*x;
    mult_sig_vis = basis_function_vis.*x_vis;
    
    basis_function2 = sin(2*pi*k*n/N); %sine wave basis function
    basis_function_vis2 = sin(2*pi*k*n_vis/(N*interp_val));
    
    mult_sig2 = basis_function2.*x;
    mult_sig_vis2 = basis_function_vis2.*x_vis;
    
    subplot(num_rows,2,(2*k+2)*2-1)
    plot(basis_function_vis)
    hold on
    plot(n*interp_val,basis_function ,'r.')
    hold off
    set(gca, 'XTick', 0:interp_val:(N)*interp_val)
    label = '';
    set(gca, 'XTickLabel', label)
    set(gca, 'Box', 'off')
    if(k > 0)
        subplot(num_rows,2,(2*k+1)*2-1)
        plot(basis_function_vis2)
        hold on
        plot(n*interp_val,basis_function2 ,'r.')
        hold off
        set(gca, 'XTick', 0:interp_val:(N)*interp_val)
        label = '';
        set(gca, 'XTickLabel', label)
        set(gca, 'Box', 'off')
    end
    
    subplot(num_rows,2,(2*k+2)*2)
    hold on
    plot(mult_sig_vis)
    plot(zeros(1, N*interp_val),'k')
    plot(n*interp_val,mult_sig ,'r.')
    hold off
    set(gca, 'XTick', 0:interp_val:(N)*interp_val)
    set(gca, 'YTick', y_range)
    set(gca, 'XTickLabel', label)
    set(gca, 'Box', 'off')
    ylim(y_range)
    xlim([0 N*interp_val])
    ylabel(num2str(round(sum(mult_sig)*100)/100))
    if(k > 0)
        subplot(num_rows,2,(2*k+1)*2)
        hold on
        plot(mult_sig_vis2)
        plot(zeros(1, N*interp_val),'k')
        plot(n*interp_val,mult_sig2 ,'r.')
        hold off
        set(gca, 'XTick', 0:interp_val:(N)*interp_val)
        set(gca, 'XTickLabel', label)
        set(gca, 'Box', 'off')
        set(gca, 'YTick', y_range)
        ylim(y_range)
        xlim([0 N*interp_val])
        ylabel(num2str(round(sum(mult_sig2)*100)/100))
    end
end

Advertisement

Filter Design Using Matlab Demo

This is the code I use in my youtube videos on filtering (ddorran channel)

%generate the noisy signal which will be filtered
x= cos(2*pi*12*[0:0.001:1.23]);
x(end) = [];
[b a] = butter(2,[0.6 0.7],'bandpass');
filtered_noise = filter(b,a,randn(1, length(x)*2));
x = (x + 0.5*filtered_noise(500:500+length(x)-1))/length(x)*2;
% This script is available at https://dadorran.wordpress.com search for
% filtering matlab demo
plot(x)
title('Noisy signal')
xlabel('Samples');
ylabel('Amplitude')

%plot magnitude spectrum of the signal
X_mags = abs(fft(x));
figure(10)
plot(X_mags)
xlabel('DFT Bins')
ylabel('Magnitude')

%plot first half of DFT (normalised frequency)
num_bins = length(X_mags);
plot([0:1/(num_bins/2 -1):1], X_mags(1:num_bins/2))
xlabel('Normalised frequency (\pi rads/sample)')
ylabel('Magnitude')

%design a second order filter using a butterworth design technique 
[b a] = butter(2, 0.3, 'low')

%Show the signal flow diagram of the filter/system
%this function available from https://dadorran.wordpress.com
create_signal_flow(b,a) 

%plot the frequency response (normalised frequency)
H = freqz(b,a, floor(num_bins/2));
hold on
plot([0:1/(num_bins/2 -1):1], abs(H),'r');

%filter the signal using the b and a coefficients obtained from
%the butter filter design function
x_filtered = filter(b,a,x);

%plot the filtered signal
figure(2)
plot(x_filtered,'r')
title('Filtered Signal - Using Second Order Butterworth')
xlabel('Samples');
ylabel('Amplitude')

%Redesign the filter using a higher order filter
[b2 a2] = butter(20, 0.3, 'low')

%Plot the magnitude spectrum and compare with lower order filter
H2 = freqz(b2,a2, floor(num_bins/2));
figure(10)
hold on
plot([0:1/(num_bins/2 -1):1], abs(H2),'g');

%filter the noisy signal and plot the result
x_filtered2 = filter(b2,a2,x);
figure(3)
plot(x_filtered2,'g')
title('Filtered Signal - Using 20th Order Butterworth')
xlabel('Samples');
ylabel('Amplitude')

%Use a band reject filter in place of a low pass filter
[b_stop a_stop] = butter(20, [0.5 0.8], 'stop');

%plot the magnitude spectrum
H_stopband = freqz(b_stop,a_stop, floor(num_bins/2));
figure(10)
hold on
plot([0:1/(num_bins/2 -1):1], abs(H_stopband),'c');

%plot filtered signal
x_filtered_stop = filter(b_stop,a_stop,x);
figure(4);
plot(x_filtered_stop,'c')
title('Filtered Signal - Using Stopband')
xlabel('Samples');
ylabel('Amplitude')

%Use matlabs built-in buttord function to get the optimum order to meet a specification
[N Wn] = buttord(0.1, 0.5, 5, 40)

%use the N and Wn values obtained above to design the filter in the usual way
[b3 a3] = butter(N, Wn, 'low');

%plot the magnitude spectrum
H3 = freqz(b3,a3, floor(num_bins/2));
figure(10);
hold on
plot([0:1/(num_bins/2 -1):1], abs(H2),'k');
figure(10)

%filter the signal and plot the ouput of the filter
x_filtered3 = filter(b3,a3,x);
figure(5);
plot(x_filtered3,'k')
title(['Filtered Signal - Using ' num2str(N) ' th Order Butterworth'])
xlabel('Samples');
ylabel('Amplitude')
% script available from https://dadorran.wordpress.com
% comparison with other filter design techniques (chebyshev and elliptical)
[b_butter a_butter] = butter(4, 0.2, 'low');
H_butter = freqz(b_butter, a_butter);

[b_cheby a_cheby] = cheby1(4, 0.5, 0.2, 'low');
H_cheby = freqz(b_cheby, a_cheby);

[b_ellip a_ellip] = ellip(4, 0.5, 40, 0.2, 'low');
H_ellip = freqz(b_ellip, a_ellip);

%plot each filter to compare 
figure(11)
norm_freq_axis = [0:1/(512 -1):1];
plot(norm_freq_axis, abs(H_butter))
hold on
plot(norm_freq_axis, abs(H_cheby),'r')
plot(norm_freq_axis, abs(H_ellip),'g')
legend('Butterworth', 'Chebyshev', 'Elliptical')
xlabel('Normalised Frequency');
ylabel('Magnitude')

%plot in dB for verification that spec is met
figure(12);
plot(norm_freq_axis, 20*log10(abs(H_butter)))
hold on
plot(norm_freq_axis, 20*log10(abs(H_cheby)),'r')
plot(norm_freq_axis, 20*log10(abs(H_ellip)),'g')
legend('Butterworth', 'Chebyshev', 'Elliptical')
xlabel('Normalised Frequency ');
ylabel('Magnitude (dB)')

create_filter_spec_plot

% Creates a plot showing a filter specification given cutoff frequencies, passband ripple and stopband attenuation.
% Optionally displays the frequency response of a filter superimposed if give b and a coefficients of the filter.
%
% Frequencies given normalised to between 0 and 1. 1 equates to Nyquist frequency or pi/rads per sample. 
%
% Example 1. Low pass filter with passband cutoff at 0.3, passband ripple of 0.6 dB, stopband at 0.8 with stopband attenuation of 20db.
% create_filter_spec_plot(0.3,0.8, 0.3 ,20)
%
% Example 2. Stop band filter with stopband edge frequencies of 0.3 and 0.5. Passband edge frequencies of 0.1 and 0.8. Passband ripple of 4dB and stopband attenuation of 50dB
% create_filter_spec_plot([0.1 0.8],[0.3 0.5], 4 ,50)
%
% Example 3. Demonstrate how to display filter frequency response with filter spec.
%  [b a] = butter(5, 1.2, 'low');
%  create_filter_spec_plot(0.1, 0.5, 0.3, 40,b,a)
function create_filter_spec_plot(PB_edge_freq,SB_edge_freq, PB_ripple ,SB_attenuation, varargin)



if nargin == 6

    b = varargin{1}

    a = varargin{2}

else

    b = [];

    a = [];

end



    



PB_edge_freq = round(PB_edge_freq*1000)/1000;

SB_edge_freq = round(SB_edge_freq*1000)/1000;

if length(PB_edge_freq) == 1 && length(SB_edge_freq) == 1

    %low pass or high pass

    if(PB_edge_freq > SB_edge_freq)

        filter_type = 'High Pass';

        PB(1,:) = [PB_edge_freq 1];

        SB(1,:) = [0 SB_edge_freq];

    else

        filter_type = 'Low Pass';

        SB(1,:) = [SB_edge_freq 1];

        PB(1,:) = [0 PB_edge_freq];

    end

else

    % band pass or stop band

    if(PB_edge_freq(1) > SB_edge_freq(1))

        filter_type = 'Band Pass';

        PB(1,:) = [PB_edge_freq(1) PB_edge_freq(2)];

        SB(1,:) = [0 SB_edge_freq(1)];

        SB(2,:) = [SB_edge_freq(2) 1];

    else

        filter_type = 'Stop Band';

        SB(1,:) = [SB_edge_freq(1) SB_edge_freq(2)];

        PB(1,:) = [0 PB_edge_freq(1)];

        PB(2,:) = [PB_edge_freq(2) 1];

    end

end





figure

hold on

num_bins = 10000;

fax = 0:1/(num_bins-1):1;





    lowest_display_attenuation = -1*(SB_attenuation*2);

[rows cols] = size(PB);

for k = 1: rows

    sig1 = ones(1, num_bins)*NaN;

    sig1(round(PB(k,1)*num_bins)+1 : round(PB(k,2)*num_bins)) = 0;

    

    sig2 = ones(1, num_bins)*NaN;

    sig2(round(PB(k,1)*num_bins)+1 : round(PB(k,2)*num_bins)) = -1*PB_ripple;



    

    sig2(round(PB(k,1)*num_bins)+1) = sig1(round(PB(k,1)*num_bins)+1);

    sig1(round(PB(k,2)*num_bins)) = sig2(round(PB(k,2)*num_bins));

    plot(fax,sig1,'k');

    plot(fax,sig2,'k');

    

    for m = 0-PB_ripple/4 : -PB_ripple/4 : -3*PB_ripple/4;

        sig3 =  ones(1, num_bins)*NaN;

        sig3(round(PB(k,1)*num_bins)+1 : round(PB(k,2)*num_bins)) = m;

        plot(fax,sig3,'k');

    end

end





[rows cols] = size(SB);



for k = 1: rows

    num_bins

    sig1 = ones(1, num_bins)*NaN;

    sig1(round(SB(k,1)*num_bins)+1 : round(SB(k,2)*num_bins)) = lowest_display_attenuation;

    

    sig2 = ones(1, num_bins)*NaN;

    sig2(round(SB(k,1)*num_bins)+1 : round(SB(k,2)*num_bins)) = -1*SB_attenuation;

    sig2(round(SB(k,1)*num_bins)+1 ) = sig1(round(SB(k,1)*num_bins)+1);

    sig1(round(SB(k,2)*num_bins)) = sig2(round(SB(k,2)*num_bins));

    plot(fax,sig1,'k');

    plot(fax,sig2,'k');

    

    sig3 = ones(1, num_bins)*NaN;

    sig3(round(SB(k,1)*num_bins)+1 : round(SB(k,2)*num_bins)) = lowest_display_attenuation: (-1*SB_attenuation-lowest_display_attenuation)/(length(round(SB(k,1)*num_bins)+1 : round(SB(k,2)*num_bins))-1) :-1*SB_attenuation;

    plot(fax,sig3,'k');

    sig4 = ones(1, num_bins)*NaN;

    sig4(round(SB(k,1)*num_bins)+1 : round(SB(k,2)*num_bins)) = fliplr(lowest_display_attenuation: (-1*SB_attenuation-lowest_display_attenuation)/(length(round(SB(k,1)*num_bins)+1 : round(SB(k,2)*num_bins))-1) :-1*SB_attenuation);

    plot(fax,sig4,'k');

end



if(length(b))

    H = freqz(b,a, num_bins);

    indices = find(H==0);

    H(indices) = 10^-16; %deal with problem of log(0)

    H_db = 20*log10(abs(H));

    %lowest_display_attenuation = min(20*log10(abs(H)))*1.2;

    plot(fax, H_db);



end



hold off

ylim([lowest_display_attenuation 0 ]);

xlim([0 1 ]);

ylabel('System Gain (db)');

xlabel('Normalised Frequency (x\pi rads/sample)')



%title([filter_type ' filter. Pass band ripple = ' num2str(PB_ripple) 'dB; Stop band attenuation = ' num2str(SB_attenuation) 'dB.' ]);