Define original function

clear
N = 2^16;
period = 1;
t = (0:N-1)*(period/N);
dt = t(2)-t(1);
fs = 1/dt;
fprintf('sample frequency %g kHz\n',fs*1e-3);
freq = [100;150;500];
A = [2, 2, 2];
nf = length(freq); fmax = 2*max(freq);
% A is 1x3, freq is 3x1 t is 1xN
x = A*cos(2*pi*freq*t);
x_in = x;
idx = 1:2:(50e-3/dt);
plot(t(idx)*1E3,x_in(idx));
xlabel('time (msec)');
ylabel('signal');

sample frequency 65.536 kHz

Fourier transform

fourier_xform(x_in,t);

inverse Fourier transform

[fx ft] = fourier_xform(x_in,t);
inv_fourier_xform(fx,ft);

maximum imaginary part 0

test round-trip difference

x_out = inv_fourier_xform(fx,ft);
diff = x_out - x;
fprintf('maximum round-trip difference: %g\n',max(diff));

maximum round-trip difference: 4.44089e-015

define lowpass filter

RC = 7.07e-4; % fc = 225.Hz
hf1 = 1./(1+j*2*pi*ft*RC);
idx = find(abs(ft)<500);
plot(ft(idx),abs(hf1(idx)));
xlabel('frequency');
ylabel('response');
hold on
r = 1/sqrt(2);
plot([-500 500], [r r],'k--',[225 225],[0.4 1],'k--');
hold off;

apply filter to input

fx1 = fx.*hf1;
inv_fourier_xform(fx1,ft);

maximum imaginary part 1.16629e-015

original signal - with noise

noi =0.5*randn(1,N);
x_noise = x+noi;
idx = 1:2:(50e-3/dt);
plot(t(idx)*1E3,x_noise(idx));
xlabel('time (msec)');
ylabel('signal');

noise spectrum

fourier_xform(noi,t);

signal and noise spectrum

fourier_xform(x_noise,t);

apply filter to noise input

[fxn ft] = fourier_xform(x_noise,t);
fx2 = fxn.*hf1;
inv_fourier_xform(fx2,ft);

maximum imaginary part 9.40238e-006

Published with MATLAB® 7.6