| name | matlab-analyze-spectrum |
| description | Analyze signal spectra in MATLAB using nonparametric and parametric estimators. Use when computing frequency content, PSD, power spectrum, spectral peaks, bandwidth, or streaming spectral analysis. Covers pspectrum, pwelch, periodogram, pmtm, pburg, pmusic, rootmusic, plomb, poctave, spectrumAnalyzer, dsp.SpectrumEstimator. TRIGGER: FFT, PSD, power spectrum, spectral analysis, frequency content, bandwidth measurement, spectral leakage, window selection, streaming spectrum, periodicity, resolve close frequencies, multitaper. DO NOT TRIGGER: filter design, spectrogram/time-frequency, audio features.
|
| license | https://www.mathworks.com/content/dam/mathworks/license/pmrl/license.md |
| metadata | {"author":"MathWorks","version":"1.0"} |
Analyze Signal Spectra in MATLAB
Compute, visualize, and extract measurements from signal spectra using Signal Processing Toolbox. Choose the right spectral estimator, avoid common FFT pitfalls, and use built-in measurement functions instead of reinventing them.
When to Use
- Computing frequency content of a signal (FFT, periodogram, PSD, power spectrum)
- Estimating power spectral density (Welch, periodogram, parametric)
- Comparing frequency content of two signals (cross-spectrum, coherence)
- Finding spectral peaks (dominant frequencies, harmonics)
- Measuring bandwidth, band power, or occupied bandwidth
- Visualizing frequency-domain data
- Streaming (real-time, frame-by-frame) spectral analysis with
spectrumAnalyzer, dsp.SpectrumEstimator, dsp.CrossSpectrumEstimator
When NOT to Use
- Designing or applying filters -- use
matlab-design-digital-filter
- Time-frequency/spectrogram analysis (including
spectrumAnalyzer spectrogram ViewType, STFT, reassignment, CWT) -- use the time-frequency analysis skill (TBD)
- Audio-specific features (MFCC, pitch, auditory models) -- use Audio Toolbox
- Control system frequency response (Bode, Nyquist) -- use Control System Toolbox
Workflow
- Batch or streaming? -- Determine if data is already in workspace (batch) or arrives frame-by-frame (streaming). See below.
- Preliminary analysis -- When user provides data, run diagnostics to inform technique selection (see below)
- Pre-process -- Check for NaN/gaps (use
plomb for non-uniform data). Apply detrend(x) only if DC/trend visibly distorts the spectrum
- Choose output type -- PSD, power spectrum, amplitude spectrum, or octave bands (see below)
- Choose estimator -- Use preliminary analysis results + Estimator Selection below
- Compute spectrum -- Use the chosen function with explicit
Fs (if available)
- Visualize -- Plot with proper axis labels (Hz, dB or linear units)
- Extract metrics -- Use
findpeaks, bandpower, obw, powerbw. See references/peak-detection.md
- Verify -- Parseval's theorem, known tone levels. See
references/power-from-spectrum.md
Batch vs Streaming
Determine the processing mode before choosing tools:
| Indicator | Mode | Tools |
|---|
| Entire signal in workspace (vector/matrix) | Batch | pwelch, pspectrum, periodogram, pburg, etc. |
| Data arrives frame-by-frame (sensor, ADC, real-time loop) | Streaming | spectrumAnalyzer, dsp.SpectrumEstimator, dsp.CrossSpectrumEstimator |
| Need live-updating spectrum display | Streaming | spectrumAnalyzer (ViewType: spectrum) |
| Need spectrum values for real-time decisions | Streaming | dsp.SpectrumEstimator (returns numeric output each frame) |
If streaming: See references/streaming-spectral-analysis.md for System object APIs, properties, and patterns. DSP System Toolbox is required.
The rest of this skill covers batch processing. For streaming, the reference file provides full guidance including the batch-to-streaming mapping table.
Preliminary Data Analysis
When the user provides signal data (variable, file, or description), run diagnostic checks before recommending a technique. See references/preliminary-analysis.md for full code.
Quick checks (run in order — stop early if decisive):
- Sampling type — Fs given → uniform. Time vector with small jitter (<5%) →
pspectrum(x,t) (resamples internally). Large jitter or truly non-uniform → plomb
- Data quality — NaN/Inf, DC offset, trend → preprocess first
- Length — <128: parametric only. 128–1024:
periodogram/pburg. >10k: pwelch preferred
- Spectral character — Run
periodogram (Hann), compute spectral flatness + peak count + strongest peak width:
- Flatness >0.3 → wideband →
pwelch (peaks are just noise fluctuations)
- Flatness <0.1, strongest peak narrow (width < 10×df) → tonal →
periodogram(...,'power') or rootmusic
- Flatness <0.1, strongest peak wide (width ≥ 10×df) → wideband-structured (chirps, shaped bands) →
pwelch + obw
- Flatness 0.1–0.3, peaks present → mixed →
pwelch + findpeaks
- AR fitness — Fit AR models with
arburg, track prediction error vs order. If error plateaus <10% at low order → pburg is strong candidate
- Tone count — If tonal: count peaks with
findpeaks, cross-check with eigenvalue drop. Known count + resolution-limited → rootmusic(x,2*nTones,Fs)
Report findings to the user before computing the spectrum:
Preliminary Analysis:
- Sampling: uniform at [Fs] Hz
- Duration: [X] s ([N] samples), resolution: [df] Hz
- Spectral character: [wideband/tonal/mixed]
- AR fit: [quality] (order [N], residual [X]%)
Recommendation: [estimator] because [justification]
Output Type Selection
| User's Goal | Output Type | Function / Option |
|---|
| Compare to a standard or another signal (length-independent) | PSD (power/Hz) | pwelch (default), periodogram |
| Read amplitude of a specific tone | Power spectrum → sqrt | pspectrum(x,Fs) (default); periodogram(x,w,N,Fs,'power') or pwelch(x,w,[],[],Fs,'power') for explicit control |
| Measure total power in a band | PSD then integrate, or use bandpower | bandpower(x,Fs,[fLow fHigh]) |
| Octave-band levels (acoustics, vibration) | Fractional-octave spectrum | poctave(x,Fs,'BandsPerOctave',3) |
| Detect frequencies in nearly-uniform data (small jitter) | Power spectrum | pspectrum(x,t) — resamples internally |
| Detect frequencies in non-uniform/gapped data | Lomb-Scargle | plomb(x,t) |
PSD vs Power Spectrum: PSD (power/Hz) is resolution-independent — broadband noise stays flat when you change parameters. Power spectrum (power/bin) is tone-friendly — peak height equals true tone power. Use PSD for noise characterization and specs (g²/Hz, dBm/Hz); use power spectrum for reading tone amplitudes. Convert: PS = PSD × RBW where RBW = enbw(win)*Fs/segLen.
Estimator Selection
Decision tree (simple → specialized):
- Just want to see what frequencies are present? →
pspectrum(x,Fs)
- Need a smooth, low-variance PSD? →
pwelch (favored for long signals — segment averaging reduces variance) or periodogram (short signal, cannot segment)
- Need accurate tone amplitude from peaks? →
periodogram(x,win,N,Fs,'power') then sqrt
- Need octave-band levels (acoustics, vibration)? →
poctave(x,Fs,'BandsPerOctave',3)
- Comparing two signals at each frequency? →
mscohere (related?) / cpsd (shared content + phase)
- Data is nearly uniform (small jitter)? →
pspectrum(x,t) (resamples internally). Truly non-uniform or large gaps? → plomb(x,t)
- Need to resolve sinusoids closer than Fs/N? →
rootmusic / pmusic (must know # of sinusoids)
- Short data, want low-variance PSD without segmenting? →
pmtm (multitaper — uses orthogonal DPSS tapers on full record; no AR assumption)
- Short data, want smooth sidelobe-free PSD? →
pburg (assumes signal is AR-like)
Quick reference table:
| Goal | Function | When to use |
|---|
| Quick look at frequency content | pspectrum | Default -- good defaults, tune Leakage (0–1) for resolution vs sidelobes |
| Single-sided amplitude spectrum | periodogram with 'power' | Need amplitude per bin |
| PSD with low variance | pwelch | Long signals — averaging segments reduces variance; preferred when signal is long enough to segment |
| PSD of short signal, low variance | pmtm | Multitaper — uses full record with orthogonal DPSS tapers; no segmentation needed |
| PSD of short signal | periodogram | Too short to segment; high variance (single window) |
| Resolve closely-spaced sinusoids | pmusic, rootmusic | Known # of sinusoids, super-resolution |
| Smooth PSD of short AR-like data | pburg | No sidelobes; NOT for resolving tones |
Default: Start with pspectrum. Drop to pwelch/periodogram for explicit control. Never use manual FFT unless user explicitly asks — redirect to built-in functions.
pspectrum vs pwelch: pspectrum outputs power spectrum (Kaiser beta=20, ~76% overlap). pwelch defaults to PSD. They match within ~0.2 dB using equivalent parameters. See references/pspectrum-equivalence.md for exact mapping. Resolution control in pspectrum: Use FrequencyResolution (Hz) to set target resolution directly, or Leakage (0–1) to control the resolution/sidelobe trade-off indirectly. TimeResolution is spectrogram-only — do NOT use it for power spectrum. Leakage parameter: pspectrum(x,Fs,Leakage=L) where L ∈ [0,1] — 0 = minimum leakage (max sidelobe suppression, widest main lobe), 1 = maximum leakage (rectangular window, best resolution), default 0.5.
Nonparametric vs Parametric: Use nonparametric (FFT-based) when signal is long or spectral shape is unknown. Use parametric when data is short, signal matches an AR model, or you need super-resolution. See references/estimator-comparison.md for detailed trade-offs.
Window Selection
| Goal | Window | Why |
|---|
| General-purpose PSD | Hann | Good balance, 50% overlap optimal |
| Detect weak tone near strong tone | Blackman-Harris or Kaiser (beta≥10) | Low sidelobes (-92 dB) |
| Resolve two close frequencies | Hann or Kaiser (beta=5) | Narrow main lobe |
| Accurate amplitude measurement | Flat Top | ~0.01 dB scalloping loss |
| Tunable trade-off | Kaiser(N,beta) | Single parameter controls curve |
Key rule: Higher sidelobe suppression costs wider main lobe. Only suppress as much as dynamic range requires. See references/spectral-windows.md for ENBW values and overlap guidelines.
Key Functions
| Function | Purpose |
|---|
pspectrum | Power spectrum with automatic defaults; Leakage parameter (0–1) controls resolution vs sidelobe suppression |
periodogram | PSD or power spectrum via single-window DFT |
pwelch | PSD via Welch's averaged periodogram |
pmtm | PSD via Thomson's multitaper method (DPSS/Slepian tapers) |
cpsd | Cross power spectral density |
mscohere | Magnitude-squared coherence (0 to 1) |
pburg / pyulear | AR PSD (Burg / Yule-Walker) |
pmusic / peig | Pseudospectrum (MUSIC / eigenvector) |
rootmusic | Frequency estimation via root-MUSIC |
findpeaks | Locate spectral peaks |
refinepeaks | Sub-bin peak frequency/amplitude refinement (use after findpeaks) |
bandpower | Power in a frequency band |
obw / powerbw | Occupied bandwidth / 3-dB power bandwidth |
meanfreq / medfreq | Mean / median frequency of spectrum |
sfdr | Spurious free dynamic range |
toi | Third-order intercept point (two-tone intermodulation) |
spectralFlatness | Spectral flatness (0=tonal, 1=white noise) |
spectralCrest | Spectral peak-to-mean ratio |
spectralKurtosis | Spectral kurtosis (transient/non-stationarity detection) |
poctave | Fractional-octave spectrum |
plomb |
Patterns
Convention reminders (apply to ALL code below):
pow2db(x) for dB — NEVER 10*log10(x)detrend(x) only when DC offset or trend is clearly an issue (e.g., large DC spike dominates spectrum)
Single-Sided Amplitude Spectrum
Fs = 1000;
N = length(x);
win = hann(N);
% periodogram with 'power' gives single-sided power per bin
% NFFT=N keeps df=Fs/N ([] defaults to nextpow2 which shifts grid)
[pxx,f] = periodogram(x,win,N,Fs,'power');
ampSpectrum = sqrt(pxx); % power → RMS amplitude
figure
plot(f,ampSpectrum)
xlabel("Frequency (Hz)")
ylabel("Amplitude (RMS)")
title("Single-Sided Amplitude Spectrum")
Key points:
- Window correction handled automatically by
periodogram
- Frequency resolution is always
Fs/N — zero-padding only interpolates
- For manual FFT: multiply by
2/N for single-sided amplitude, do NOT double DC
PSD with Welch's Method
Fs = 44100;
% x = detrend(x); % uncomment if DC/trend distorts low-frequency estimates
segmentLength = round(Fs * 0.05); % 50 ms → df = 20 Hz
win = hann(segmentLength);
[pxx,f] = pwelch(x,win,round(0.5*segmentLength),[],Fs);
figure
plot(f,pow2db(pxx))
xlabel("Frequency (Hz)")
ylabel("PSD (dB/Hz)")
title("Power Spectral Density (Welch)")
grid on
Key points:
- Segment length sets resolution:
df = Fs / segmentLength
- More segments = lower variance but coarser resolution (trade-off)
- 50% overlap with Hann is optimal
Spectral Peak Detection
See references/peak-detection.md for noise-floor-relative thresholds and harmonic detection.
[pxx,f] = pwelch(x,hann(1024),[],[],Fs);
pxxDB = pow2db(pxx);
[peaks,locs] = findpeaks(pxxDB,f, ...
MinPeakHeight=-20, ...
MinPeakDistance=50, ...
MinPeakProminence=10);
figure
plot(f,pxxDB)
hold on
plot(locs,peaks,"rv",MarkerFaceColor="r")
hold off
xlabel("Frequency (Hz)")
ylabel("Power/Frequency (dB/Hz)")
legend("PSD","Peaks")
Band Power and Bandwidth
See references/spectral-measurements.md for full details including PSD input, power spectrum + rbw input, and spectral entropy.
Fs = 1000;
pTotal = bandpower(x,Fs,[0 Fs/2]);
pBand = bandpower(x,Fs,[50 150]);
fracPower = pBand / pTotal;
bw = obw(x,Fs); % 99% occupied bandwidth
pbw = powerbw(x,Fs); % 3-dB bandwidth (narrowband/resonance signals only)
Fractional-Octave Spectrum (poctave)
See references/poctave-api.md for full syntax, common mistakes, and octave smoothing pattern.
Fs = 44100;
% 1/3-octave with A-weighting (environmental noise)
[p,cf] = poctave(x,Fs,'BandsPerOctave',3,'Weighting','A','FrequencyLimits',[20 20000]);
% From existing PSD (no recomputation needed)
[pxx,f] = pwelch(x,[],[],[],Fs);
[p,cf] = poctave(pxx,Fs,f,'BandsPerOctave',3);
Key points:
- Use
'BandsPerOctave',3 — do NOT use '1/3 octave' (invalid syntax)
- Only valid types:
'power' (default), 'spectrogram' — NO 'psd' type
- Can accept pre-computed PSD:
poctave(pxx,Fs,f) avoids reprocessing
- Output is always power — use
pow2db(p) for dB display
Cross-Spectral Comparison
Fs = 1000;
segLen = 1024;
win = hann(segLen);
[pxy,f] = cpsd(x,y,win,segLen/2,[],Fs);
[cxy,f] = mscohere(x,y,win,segLen/2,[],Fs);
figure
tiledlayout(2,1)
nexttile
plot(f,pow2db(abs(pxy)))
ylabel("|CPSD| (dB)"); xlabel("Frequency (Hz)")
nexttile
plot(f,cxy)
ylabel("Coherence"); xlabel("Frequency (Hz)"); ylim([0 1])
Key points:
cpsd is complex — use abs for magnitude, angle for phase lag- Coherence ≈ 1: linearly related; ≈ 0: unrelated at that frequency
- Use same
win/overlap as pwelch for consistent resolution
Multitaper PSD (Thomson's Method)
Fs = 1000;
nw = 4; % time-halfbandwidth product (default); resolution bandwidth = 2*nw*Fs/N
[pxx,f] = pmtm(x,nw,[],Fs);
plot(f,pow2db(pxx))
xlabel("Frequency (Hz)"); ylabel("PSD (dB/Hz)")
Key points:
- Uses 2*nw-1 orthogonal DPSS (Slepian) tapers on the full record — no segmentation
- Variance reduction without sacrificing frequency resolution (unlike Welch)
- Resolution bandwidth =
2*nw*Fs/N — increase nw for smoother PSD (wider bandwidth), decrease for finer resolution
- Ideal for short records where Welch can't segment, or when you need maximum spectral information from limited data
- Output is true PSD — directly comparable to
pwelch
- Confidence intervals available:
[pxx,f,pxxc] = pmtm(x,nw,[],Fs,ConfidenceLevel=0.95)
- When to prefer over
pwelch: short data, can't afford to segment, or need optimal bias-variance without AR assumptions
- When to prefer over
pburg: no prior knowledge that signal is AR-like; want a nonparametric estimate
Parametric PSD (Burg Method)
Fs = 1000;
modelOrder = 16; % 2*numPeaks to 4*numPeaks
[pxx,f] = pburg(x,modelOrder,[],Fs);
plot(f,pow2db(pxx))
xlabel("Frequency (Hz)"); ylabel("PSD (dB/Hz)")
Key points:
- Order too low = over-smoothed; too high = spurious peaks
- Burg preferred over Yule-Walker (less biased for short data)
- Output is true PSD — comparable to
pwelch
Frequency Estimation with MUSIC
Fs = 1000;
nSinusoids = 2;
p = 2*nSinusoids; % subspace dim: 2 per real sinusoid (pos + neg freq)
[pseudoSpectrum,f] = pmusic(x,p,1024,Fs);
freqEstimates = rootmusic(x,p,Fs); % more accurate than reading peaks
Key points:
- MUSIC output is a pseudospectrum, NOT a PSD — peak heights ≠ power
p = 2 × number of real sinusoids — using p = nSinusoids gives wrong resultsrootmusic gives more accurate frequencies than reading pmusic peaks- Works best for sinusoids in white noise; short data OK if SNR adequate
Manual FFT-Based PSD (When User Insists)
Only when user explicitly needs manual FFT. See references/manual-fft-psd.md for windowed, power spectrum, and complex-valued variants.
Formula: psdx = (1/(Fs*N)) * abs(fft(x)(1:N/2+1)).^2, double interior bins only (not DC/Nyquist). Windowed: normalize by 1/(Fs*sum(win.^2)).
Conventions
- Use
pspectrum as default — redirect manual FFT users to built-in functions
- Detrend when warranted — apply
detrend(x) when DC/trend visibly dominates the spectrum or distorts low-frequency estimates; do not blindly detrend every signal
- Label axes with units — "Frequency (Hz)", "PSD (dB/Hz)", "Amplitude (V)"
pow2db/db2pow for dB — never 10*log10(...) manually. See references/units-and-db.md for domain-specific reference levels- Frequency resolution is
Fs/N — zero-padding only interpolates
- Use built-in measurement functions —
bandpower, obw, powerbw, findpeaks
- Favor convenience plots — calling
pwelch(x,[],[],[],Fs), pspectrum(x,Fs), pburg(x,order,[],Fs), etc. without output arguments generates a labeled plot automatically; prefer this over manual plotting when a quick visualization is needed
- Write a script — when performing multi-step analysis, write the complete workflow as a
.m script and run it with run_matlab_file, rather than executing code piecemeal across multiple evaluate_matlab_code calls. This produces a reviewable, reproducible artifact
Common Mistakes
| Mistake | Correct Approach |
|---|
fft(x) without window | Apply hann(N) or similar |
abs(fft(x)).^2 as "PSD" | Use periodogram/pwelch; if manual, scale by 1/(Fs*N) and double interior only |
| Confusing PSD with power spectrum | PSD = power/Hz; PS = power/bin |
| Zero-padded FFT = higher resolution | Zero-padding interpolates, doesn't resolve |
10*log10(pxx) | Use pow2db |
| Manual integration for band power | Use bandpower(x,Fs,[fLow fHigh]) |
Comparing pspectrum to pwelch directly | Different units (PS vs PSD). See references/pspectrum-equivalence.md |
Using pentropy | Deprecated — use spectralEntropy instead |
pspectrum(...,'TimeResolution',val) for power spectrum | TimeResolution is spectrogram-only. Use 'FrequencyResolution',val or 'Leakage',val |
obw(x,Fs,[],'Percentpower',pct) | obw uses positional args, not name-value. Correct: obw(x,Fs,[fLow fHigh],pct) where pct is the 4th positional arg |
| Filter-bank tone power reads 1–3 dB low | Tone is between bins — filter-bank has narrow pass-bands. Fix: coherent frame length (frameLen = Fs/df) so tones land on bin centers |
snr(x,Fs) unstable for multi-tone with harmonic relationships | snr classifies harmonics specially — 200 Hz = 2nd harmonic of 100 Hz causes inconsistent results. Use per-tone SNR from spectrum instead |
Troubleshooting
See references/troubleshooting.md for detailed explanations and code examples.
- "Two tones but only one peak" — Resolution limit (
Fs/N) or sidelobe masking. Fix: more data, lower-sidelobe window (blackmanharris, kaiser(N,15)), or reassigned periodogram ([pxx,f,rpxx,fc] = periodogram(x,win,N,Fs,'reassigned') — 4-output syntax returns reassigned PSD and center-of-energy frequencies for sharper peak localization).
- "Tone power is wrong" — Scalloping loss (tone between bins). Fix:
flattopwin(N) for amplitude accuracy. Worst case: Hann loses 1.4 dB, Flat Top loses 0.01 dB.
Verification
Verify with Parseval's theorem — total power from spectrum should match time-domain. See references/power-from-spectrum.md for formulas.
pTime = sum(x.^2)/length(x);
pFreq = sum(pxx) * (f(2)-f(1)); % integrate PSD
assert(abs(pTime - pFreq)/pTime < 0.02,"Parseval check failed")
Copyright 2026 The MathWorks, Inc.
