diff --git a/adc_eval/spectrum.py b/adc_eval/spectrum.py
index 3150961..66662ae 100644
--- a/adc_eval/spectrum.py
+++ b/adc_eval/spectrum.py
@@ -90,18 +90,25 @@ def find_harmonics(spectrum, freq, nfft, bin_sig, psig, harms=5, leak=20):
             bin_harm = int(harm - (zone - 1) * nfft / 2)
 
         # Make sure we pick the max bin where power is maximized; due to spectral leakage
-        bin_harm_max = bin_harm
+        # if bin_harm == nfft/2, set to bin of 0
+        if bin_harm == nfft / 2:
+            bin_harm = 0
+        pwr_max = spectrum[bin_harm]
         for i in range(bin_harm - leak, bin_harm + leak + 1):
-            if spectrum[i] > spectrum[bin_harm_max]:
-                bin_harm_max = i
-
-        bin_harm = bin_harm_max
+            try:
+                pwr = spectrum[i]
+                if pwr > pwr_max:
+                    bin_harm = i
+                    pwr_max = pwr
+            except IndexError:
+                # bin + leakage out of bounds, so stop looking
+                break
 
         harm_stats["harm"][harm_index]["bin"] = bin_harm
-        harm_stats["harm"][harm_index]["power"] = spectrum[bin_harm]
+        harm_stats["harm"][harm_index]["power"] = pwr
         harm_stats["harm"][harm_index]["freq"] = round(freq[bin_harm] / 1e6, 1)
-        harm_stats["harm"][harm_index]["dBc"] = dBW(spectrum[bin_harm] / psig)
-        harm_stats["harm"][harm_index]["dB"] = dBW(spectrum[bin_harm])
+        harm_stats["harm"][harm_index]["dBc"] = dBW(pwr / psig)
+        harm_stats["harm"][harm_index]["dB"] = dBW(pwr)
 
         harm_index = harm_index + 1
 
diff --git a/tests/test_spectrum_plotting.py b/tests/test_spectrum_plotting.py
index 5672366..d282609 100644
--- a/tests/test_spectrum_plotting.py
+++ b/tests/test_spectrum_plotting.py
@@ -70,6 +70,30 @@ def test_find_harmonics_with_leakage(self):
                     msg=msg_txt,
                 )
 
+    def test_find_harmonics_with_leakage_outside_bounds(self):
+        """Test find harmonics with leakage bins exceeding array bounds."""
+        self.bin = self.nfft / 4 - 0.5
+        (freq, pwr) = self.gen_spectrum(5)
+        leakage_bins = 2
+        stats = spectrum.find_harmonics(
+            pwr, freq, self.nfft, self.bin, self.arms, harms=2, leak=leakage_bins
+        )
+        self.assertTrue(self.nfft / 2 - 3 <= stats["harm"][2]["bin"], self.nfft / 2 - 1)
+
+    def test_find_harmonics_on_fft_bound(self):
+        """Test find harmonics with harmonics landing at nfft/2."""
+        self.nfft = 2**12
+        self.bin = self.nfft / 8
+        (freq, pwr) = self.gen_spectrum(10)
+        leakage_bins = 0
+        stats = spectrum.find_harmonics(
+            pwr, freq, self.nfft, self.bin, self.arms, harms=5, leak=leakage_bins
+        )
+        self.assertEqual(stats["harm"][2]["bin"], 2 * self.bin)
+        self.assertEqual(stats["harm"][3]["bin"], 3 * self.bin)
+        self.assertEqual(stats["harm"][4]["bin"], 0)
+        self.assertEqual(stats["harm"][5]["bin"], self.nfft - 5 * self.bin)
+
     def test_plot_string(self):
         """Test proper return of plotting string."""
         self.bin = 13