Fourier transform and LTI filter and frequency response in Matlab

Question

I'm new to Matlab for LTI signal processing and wondering if anyone can help with something that I'm sure is meant to be basic. I've spent hours and hours researching and obtaining background information and still cannot obtain a clear path to tackle these problems. So far, from scratch, I have generated a signal required and managed to use the fft function to produce the signal's DFT:

function x = fourier_rikki(A,t,O)

Fs = 1000;
t = 0:(1/Fs):1;

A = [0.5,0,0.5];
N = (length(A) - 1)/2;
x = zeros(size(t));
f1 = 85;
O1 = 2*pi*f1;

for k = 1:length(A)
x1 = x + A(k)*exp(1i*O1*t*(k-N-1));
end

f2 = 150;
O2 = 2*pi*f2;

for k = 1:length(A);
x2 = x + A(k)*exp(1i*O2*t*(k-N-1));
end

f3 = 330;
O3 = 2*pi*f3;

for k = 1:length(A);
x3 = x + A(k)*exp(1i*O3*t*(k-N-1));
end

signal = x1 + x2 + x3;

figure(1);
subplot(3,1,1);
plot(t, signal);
title('Signal x(t) in the Time Domain');
xlabel('Time (Seconds)');
ylabel('x(t)');

X = fft(signal);  %DFT of the signal

subplot(3,1,2);
plot(t, X);
title('Power Spectrum of Discrete Fourier Transform of x(t)');
xlabel('Time (Seconds)');
ylabel('Power');

f = linspace(0, 1000, length(X)); %?

subplot(3,1,3);
plot(f, abs(X));  %Only want the positive values
title('Spectral Frequency');
xlabel('Frequency (Hz)'); ylabel('Power');

end

At this stage, I'm assuming this is correct for:

"Generate a signal with frequencies 85,150,330Hz using a sampling frequency of 1000Hz - plot 1seconds worth of the signal and its Discrete Fourier Transform."

The next step is to "Find the frequency response of an LTI system that filters out the higher and lower frequencies using the Fourier Transform". I'm stuck trying to create an LTI system that does that! I have to be left with the 150Hz signal, and I'm guessing I perform the filtering on the FFT, perhaps using conv.

My course is not a programming course - we are not assessed on our programming skills and I have minimal Matlab experience - basically we have been left to our own devices to struggle through, so any help would be greatly appreciated! I am sifting through tonnes of different examples and searching Matlab functions using 'help' etc, but since each one is different and does not have a break down of the variables used, explaining why certain parameters/values are chosen etc. it is just adding to the confusion.

Among many (many) others I have looked at: http://www.ee.columbia.edu/~ronw/adst-spring2010/lectures/matlab/lecture1.html http://gribblelab.org/scicomp/09_Signals_and_sampling.html section 10.4 especially. As well as Matlab Geeks examples and Mathworks Matlab function explanations. I guess the worst that can happen is that nobody answers and I continue burning my eyeballs out until I manage to come up with something :) Thanks in advance.

I found this bandpass filter code as a Mathworks example, which is exactly what needs to be applied to my fft signal, but I don't understand the attenuation values Ast or the amount of ripple Ap.

n = 0:159;
x = cos(pi/8*n)+cos(pi/2*n)+sin(3*pi/4*n);

d = fdesign.bandpass('Fst1,Fp1,Fp2,Fst2,Ast1,Ap,Ast2',1/4,3/8,5/8,6/8,60,1,60);
Hd = design(d,'equiripple');

y = filter(Hd,x);
freq = 0:(2*pi)/length(x):pi;
xdft = fft(x);
ydft = fft(y);

plot(freq,abs(xdft(1:length(x)/2+1)));
hold on;
plot(freq,abs(ydft(1:length(x)/2+1)),'r','linewidth',2);
legend('Original Signal','Bandpass Signal');

Answer 1

Here is something you can use as a reference. I think I got the gist of what you were trying to do. Let me know if you have any questions.

clear all
close all

Fs = 1000;
t = 0:(1/Fs):1;
N = length(t);

% 85, 150, and 330 Hz converted to radian frequency
w1 = 2*pi*85;
w2 = 2*pi*150;
w3 = 2*pi*330;

% amplitudes
a1 = 1;
a2 = 1.5;
a3 = .75;

% construct time-domain signals
x1 = a1*cos(w1*t);
x2 = a2*cos(w2*t);
x3 = a3*cos(w3*t);

% superposition of 85, 150, and 330 Hz component signals
x = x1 + x2 + x3; 

figure
plot(t(1:100), x(1:100));
title('unfiltered time-domain signal, amplitude vs. time');
ylabel('amplitude');
xlabel('time (seconds)');

% compute discrete Fourier transform of time-domain signal
X = fft(x); 
Xmag = 20*log10(abs(X)); % magnitude spectrum
Xphase = 180*unwrap(angle(X))./pi; % phase spectrum (degrees)
w = 2*pi*(0:N-1)./N; % normalized radian frequency
f = w./(2*pi)*Fs; % radian frequency to Hz
k = 1:N; % bin indices 

% plot magnitude spectrum
figure
plot(f, Xmag)
title('frequency-domain signal, magnitude vs. frequency');
xlabel('frequency (Hz)');
ylabel('magnitude (dB)');

% frequency vector of the filter. attenuates undesired frequency components
% and keeps desired components. 
H = 1e-3*ones(1, length(k));
H(97:223) = 1;
H((end-223):(end-97)) = 1;

% plot magnitude spectrum of signal and filter
figure
plot(k, Xmag)
hold on
plot(k, 20*log10(H), 'r')
title('frequency-domain signal (blue) and filter (red), magnitude vs. bin index');
xlabel('bin index');
ylabel('magnitude (dB)');

% filtering in frequency domain is just multiplication
Y = X.*H;

% plot magnitude spectrum of filtered signal
figure
plot(f, 20*log10(abs(Y)))
title('filtered frequency-domain signal, magnitude vs. frequency');
xlabel('frequency (Hz)');
ylabel('magnitude (dB)');

% use inverse discrete Fourier transform to obtain the filtered time-domain
% signal. This signal is complex due to imperfect symmetry in the
% frequency-domain, however the imaginary components are nearly zero.
y = ifft(Y);

% plot overlay of filtered signal and desired signal
figure
plot(t(1:100), x(1:100), 'r')
hold on
plot(t(1:100), x2(1:100), 'linewidth', 2)
plot(t(1:100), real(y(1:100)), 'g')
title('input signal (red), desired signal (blue), signal extracted via filtering (green)');
ylabel('amplitude');
xlabel('time (seconds)');

Here is the end result...