diff --git a/virgo/virgo.py b/virgo/virgo.py
index 459f860..1f55797 100644
--- a/virgo/virgo.py
+++ b/virgo/virgo.py
@@ -567,7 +567,7 @@ def observe(obs_parameters, spectrometer='wola', obs_file='observation.dat', sta
 t_sample='''+str(t_sample)+'''
 duration='''+str(duration))
 
-def plot(obs_parameters='', n=0, m=0, f_rest=0, slope_correction=False, dB=False, rfi=[0,0], xlim=[0,0], ylim=[0,0], dm=0,
+def plot(obs_parameters='', n=0, m=0, f_rest=0, slope_correction=False, dB=False, rfi=[], xlim=[0,0], ylim=[0,0], dm=0,
 	 obs_file='observation.dat', cal_file='', waterfall_fits='', spectra_csv='', power_csv='', plot_file='plot.png'):
 	'''
 	Process, analyze and plot data.
@@ -684,14 +684,18 @@ def best_fit(power):
 	waterfall = waterfall[3:, :]
 
 	# Mask RFI-contaminated channels
-	if rfi != [0,0]:
-		# Frequency to channel transformation
-		rfi_lo = channels*(rfi[0] - (frequency - bandwidth/2))/bandwidth
-		rfi_hi = channels*(rfi[1] - (frequency - bandwidth/2))/bandwidth
+	if rfi != []:
+		
+		for j in range(len(rfi)):
+			
+			# Frequency to channel transformation
+			current_rfi = rfi[j]
+			rfi_lo = channels*(current_rfi[0] - (frequency - bandwidth/2))/bandwidth
+			rfi_hi = channels*(current_rfi[1] - (frequency - bandwidth/2))/bandwidth
 
-		# Blank channels
-		for i in range(int(rfi_lo), int(rfi_hi)):
-			waterfall[:, i] = np.nan
+			# Blank channels
+			for i in range(int(rfi_lo), int(rfi_hi)):
+				waterfall[:, i] = np.nan
 
 	if cal_file != '':
 		waterfall_cal = offset*np.fromfile(cal_file, dtype='float32').reshape(-1, channels)/bins
@@ -700,11 +704,19 @@ def best_fit(power):
 		waterfall_cal = waterfall_cal[3:, :]
 
 		# Mask RFI-contaminated channels
-		if rfi != [0,0]:
-			# Blank channels
-			for i in range(int(rfi_lo), int(rfi_hi)):
-				waterfall_cal[:, i] = np.nan
-
+		if rfi != []:
+			
+			for j in range(len(rfi)):
+				
+				# Frequency to channel transformation
+				current_rfi = rfi[j]
+				rfi_lo = channels*(current_rfi[0] - (frequency - bandwidth/2))/bandwidth
+				rfi_hi = channels*(current_rfi[1] - (frequency - bandwidth/2))/bandwidth
+			
+				# Blank channels
+				for i in range(int(rfi_lo), int(rfi_hi)):
+					waterfall_cal[:, i] = np.nan
+					
 	# Compute average spectra
 	with warnings.catch_warnings():
 		warnings.filterwarnings(action='ignore', message='Mean of empty slice')