Matlab Analyze Spectrum

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:

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

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:

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

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

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

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.