a = imread('Image.bmp'); % read the bmp file

figure(1)
image(a);     % display the colored picture

figure(2)
ad = sum(a,3); % add the three components of the RGB picture
colormap(gray)
imagesc(ad);     % display gray scale picture picture

%%% choose one line for analysis
nn = 300;
%%nn = 100;
inline = ad(nn,:);

dt = 5/10000;   % set the sampling interval, frequency 2000 Hz 
np = length(inline);
t = dt*[0:np-1];
figure(3)
plot(t,inline)
grid on
xlabel('time [sec]'); ylabel('signal');

T = np*dt;
freq = 1/T*[0:np-1];

c = fft(inline);
amp = abs(c)/np;
figure(4)
plot(freq,amp)
maxAmp = max(amp(10:np-10));
axis([0, max(freq) ,0, maxAmp])
grid on
xlabel('frequency [Hz]'); ylabel('amplitude');

maxFreq = 300;
%%nMax = find(freq>maxFreq)(1);  %% works in Octave only
bigFreq = find(freq>maxFreq);
nMax =  bigFreq(1);

figure(4)
plot(freq(1:nMax),amp(1:nMax))
axis([0, freq(nMax),0,maxAmp])
grid on
xlabel('frequency [Hz]'); ylabel('amplitude');

bigAmps = find(amp>0.8*maxAmp)

%%%% remove the component with frequency 120 Hz
n1 = 55;  c(n1) = 0; c(np+2-n1) = 0;  % central frequency
          c(n1-1) = 0; c(np+2-n1+1) = 0;  % one smaller
          c(n1+1) = 0; c(np+2-n1-1) = 0;  % one larger
signal2 = real(ifft(c));              % reconstruct signal

%%%% remove the component with frequency 2*120 Hz
n2 = 109;  c(n2) = 0; c(np+2-n2) = 0;  % central frequency

figure(3)
plot(t,inline,t,signal2)
grid on
xlabel('time [sec]'); ylabel('signal');

%%%% run the fft for all rows simultaneously
cM = fft(ad,[],2);
cM(:,n1)   = 0; cM(:,np+2-n1)= 0;
cM(:,n1-1) = 0; cM(:,np+2-n1+1)= 0;
cM(:,n1+1) = 0; cM(:,np+2-n1-1)= 0;
cM(:,n2)   = 0; cM(:,np+2-n2)= 0;
adNew = real(ifft(cM, [] ,2));

figure(5)
colormap(gray)
imagesc(adNew)

figure(6)
colormap(1-gray);
cplot = cM(:,1:200); 
image(abs(cplot)/50)
