Code covered by the BSD License  

Highlights from
basic signal processing functions

basic signal processing functions

by

 

h.fft(x,N,W) h.ifft(x,N) h.conv(x,y) h.freqz(b,a,N,W) h.interp(x,N) h.decimate(x,N)

mydsp
function h = mydsp
% implement some basic DSP related matlab built-in functions
% h.fft(x,N,W)
% h.ifft(x,N)
% h.conv(x,y)
% h.freqz(b,a,N,W)
% h.interp(x,N)
% h.decimate(x,N)
%
% by: yuanye wang
% Feb 2013

h.fft = @myfft;
h.ifft = @myifft;
h.conv = @myconv;
h.freqz = @myfreqz;
h.interp = @myinterp;
h.decimate = @mydecimate;
end

function [y,W] = myfft(x,N,W)
% X(n) = sum_1^M( x(m)*exp(-j*2*pi*(n-1)*(m-1)/N) )
if nargin == 1, N = length(x); W = 2*pi*(0:N-1)/N; end
if nargin == 2, W = 2*pi*(0:N-1)/N; end
if (length(x)>N), x = x(1:N); fprintf('input too long, truncated to length N;\n'); end
m = 0:length(x)-1;
y = x * exp(-1j*m'*W);
end

function y = myifft(x,N)
% X(m) = 1/N*sum_1^N( x(n)*exp(j*2*pi*(n-1)*(m-1)/N) )
if (length(x)>N), x = x(1:N); fprintf('input too long, truncated to length N;\n'); end
if nargin == 1, N = length(x); end
n = 0:length(x)-1;
W = 2*pi*(0:N-1)/N;
y = 1/N * x * exp(1j*n'*W);
end

function z = myconv(x,y)
% z[n] = x[n] conv y[n] = sum(x[k]*y[n-k])
N = length(x)+length(y)-1;
z = zeros(1,N);
for n = 1:N
    y_idx = n - (1:length(x)) + 1;
    y_idx(y_idx<=0) = []; y_idx(y_idx>length(y)) = [];
    z(n) = x(n-y_idx+1) * transpose(y(y_idx));
end
end

function [H,W] = myfreqz(b,a,N,W)
% H(jw) = (b0+b1*exp(-j2*pi/N*)+...)/(a0+a1*exp(-jw)+...)
if nargin == 2, N = 512; W = 2*pi*(0:N-1)/N; end
if nargin == 3, W = 2*pi*(0:N-1)/N; end
B = b * exp(-1j*(0:length(b)-1)'*W);
A = a * exp(-1j*(0:length(a)-1)'*W);
H = B./A;
if nargout == 0
    figure;
    subplot(2,1,1);
    plot(W,20*log10(abs(H))); hold on;
    xlabel('Frequency'); ylabel('Magnitude [dB]');
    title('freqz function')
    subplot(2,1,2);
    plot(W,angle(H)/pi*180); hold on;
    xlabel('Frequency'); ylabel('Phase [degree]');
end
end

function y = myinterp(x,N)
l = 4; % same as Matlab build-in func
nTap = l * N * 2;
h = fir1(nTap,1/N,'low');
delay = nTap/2;

len_x = length(x);
x = [x zeros(1,delay)];
x = [x; zeros(N-1,length(x))];
x = reshape(x(:),1,[]);
y = N*filter(h,1,x);
y = y(delay+(1:len_x*N));
end

function y = mydecimate(x,N)
l = 4; % same as Matlab build-in func
nTap = l * N * 2;
h = fir1(nTap,1/N,'low');
delay = nTap/2;

len_x = length(x);
x = [x zeros(1,delay)];
y = filter(h,1,x);
y = y(delay+(1:len_x));
y = y(1:N:end);
end

Contact us