matlab-analyze-spectrum

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

1. Batch or streaming? -- Determine if data is already in workspace (batch) or arrives frame-by-frame (streaming). See below.
2. Preliminary analysis -- When user provides data, run diagnostics to inform technique selection (see below)
3. Pre-process -- Check for NaN/gaps (use plomb for non-uniform data). Apply detrend(x) only if DC/trend visibly distorts the spectrum
4. Choose output type -- PSD, power spectrum, amplitude spectrum, or octave bands (see below)
5. Choose estimator -- Use preliminary analysis results + Estimator Selection below
6. Compute spectrum -- Use the chosen function with explicit Fs (if available)
7. Visualize -- Plot with proper axis labels (Hz, dB or linear units)
8. Extract metrics -- Use findpeaks, bandpower, obw, powerbw. See references/peak-detection.md
9. 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):

1. Sampling type — Fs given → uniform. Time vector with small jitter (<5%) → pspectrum(x,t) (resamples internally). Large jitter or truly non-uniform → plomb
2. Data quality — NaN/Inf, DC offset, trend → preprocess first
3. Length — <128: parametric only. 128–1024: periodogram/pburg. >10k: pwelch preferred
4. 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
5. AR fitness — Fit AR models with arburg, track prediction error vs order. If error plateaus <10% at low order → pburg is strong candidate
6. 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):

1. Just want to see what frequencies are present? → pspectrum(x,Fs)
2. Need a smooth, low-variance PSD? → pwelch (favored for long signals — segment averaging reduces variance) or periodogram (short signal, cannot segment)
3. Need accurate tone amplitude from peaks? → periodogram(x,win,N,Fs,'power') then sqrt
4. Need octave-band levels (acoustics, vibration)? → poctave(x,Fs,'BandsPerOctave',3)
5. Comparing two signals at each frequency? → mscohere (related?) / cpsd (shared content + phase)
6. Data is nearly uniform (small jitter)? → pspectrum(x,t) (resamples internally). Truly non-uniform or large gaps? → plomb(x,t)
7. Need to resolve sinusoids closer than Fs/N? → rootmusic / pmusic (must know # of sinusoids)
8. Short data, want low-variance PSD without segmenting? → pmtm (multitaper — uses orthogonal DPSS tapers on full record; no AR assumption)
9. 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 results
• rootmusic 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.

---