fixes filter on MPT-transition. Add a high quality filter wir anti-aliasing

- Implement cross-correlation-based phase inversion detection
- Add --force-invert-phase flag for manual override and testing
- Add --no-phase-correction flag to disable automatic detection
- Update README with comprehensive documentation
- Improve phase detection sensitivity and add detailed logging
- Ensure consistent IR polarity for easier mixing of multiple IRs

CHANGELOG.md

@@ -1,5 +1,49 @@
 # Changelog
 
+## [v1.2.1] - 2024-12-20
+### Added
+- **High-Quality Filter Implementation**: Direct Form II Transposed filter structure for improved numerical stability
+- **Filter Warm-Up**: Proper filter initialization to avoid transients, especially important for MPT transforms
+- **Anti-Aliasing in Resampling**: Automatic anti-aliasing filter when downsampling to prevent aliasing artifacts
+
+### Fixed
+- **MPT Filter Application**: MPT IRs now correctly use the filtered recorded sweep (same as RAW IRs)
+- **Filter Response at Lower Sample Rates**: Fixed filter response accuracy when resampling to 48kHz and other lower sample rates
+- **Resampling Artifacts**: Added proper anti-aliasing to prevent frequency aliasing when downsampling
+
+### Changed
+- **Filter Quality**: Improved filter implementation with high-precision coefficient calculation and proper state management
+- **Filter Artifacts**: Reduced artifacts in steep filters (48+ dB/octave) through better numerical stability
+
+### Technical Improvements
+- **Direct Form II Transposed**: More numerically stable filter structure for cascaded filters
+- **Filter Warm-Up**: Pre-initialization of filter state to eliminate transients
+- **Anti-Aliasing**: Steep low-pass filter (48 dB/octave) applied before downsampling to prevent aliasing
+
+## [v1.2.0] - 2024-12-20
+### Added
+- **Cabpack Generation**: Complete cabpack creation with IRs in 9 different formats organized in structured directory trees
+- **Default Mode**: Run without command-line options for automatic cabpack generation from current directory
+- **Executable Directory Detection**: Automatic detection of binary location for double-click support in Finder/Explorer
+- **Automatic File Organization**: Support for `selection.txt` and `ambient.txt` files to automatically organize IRs
+- **Mixes Folder Support**: Automatic conversion of ready-to-use IRs from `mixes` folder to all cabpack formats
+- **IR Filename in Plots**: IR filenames now displayed in plot titles for better identification
+- **Cabpack Defaults**: Automatic defaults for cabpack mode (fade-ms: 10, lowcut: 40Hz)
+- **Plots for Mixes**: Automatic plot generation for mix IRs in `plots/mixes` folder
+- **Logo in Plots**: Valhallir logo automatically displayed in upper-left corner of all generated plots (embedded in binary)
+
+### Changed
+- **Optional Command-Line Flags**: All main flags (`--sweep`, `--recorded`, `--output`) are now optional with intelligent defaults
+- **Default Sweep File**: Changed default sweep filename to `sweep.wav` (simpler naming)
+- **Output Directory**: IRs folder created one directory level up from recorded directory by default
+- **Directory Structure**: Added `mixes` subfolder to each format folder (only created if `mixes` folder exists)
+- **Conditional Folder Creation**: `selection`, `ambient mics`, and `mixes` folders are only created if the corresponding files/folders exist in the recorded directory
+
+### Technical Improvements
+- **Smart Path Detection**: Uses executable directory when double-clicked, current directory when run from command line
+- **Enhanced Logging**: Better progress information and path logging for debugging
+- **Improved Error Handling**: Graceful handling of missing optional files (selection.txt, ambient.txt, mixes folder)
+
 ## [v1.1.0] - 2024-12-19
 ### Added
 - **IR Visualization**: New `--plot-ir` flag to generate frequency response and waveform plots

README.md

@@ -10,10 +10,14 @@ A CLI tool for processing WAV files to generate impulse responses (IR) from swee
 - **Optional low-cut and high-cut filtering:** Apply Butterworth filters to the recorded sweep before IR extraction (--lowcut, --highcut, --cut-slope)
 - **Automatic fade-out:** Linear fade-out at the end of the IR to avoid clicks (default 5 ms, configurable with --fade-ms)
 - **IR Visualization:** Generate frequency response and waveform plots with `--plot-ir`
+- **Automatic Phase Correction:** Detects and corrects phase-inverted recorded sweeps for consistent IR polarity
+- **Manual Phase Inversion:** Use `--force-invert-phase` to explicitly invert the recorded sweep for testing or manual override
 - **96kHz 24-bit WAV file support** for high-quality audio processing
 - **Multiple output formats** with configurable sample rates and bit depths
 - **Minimum Phase Transform (MPT)** option for reduced latency IRs
 - **Automatic silence trimming** and normalization
+- **Batch processing:** Process entire directories of recorded files automatically
+- **Cabpack generation:** Create complete cabpack structures with IRs in multiple formats organized in a directory tree
 - **Modular design** with separate packages for WAV I/O, convolution, and visualization
 - **Robust error handling** and validation
 
@@ -30,6 +34,24 @@ go build -o valhallir-deconvolver
 
 ## Usage
 
+### Default Mode (No Options)
+
+If you run the tool without any command-line options, it automatically creates a cabpack:
+
+```sh
+./valhallir-deconvolver
+```
+
+This will:
+- Use the current directory as the recorded directory
+- Look for `sweep.wav` in the current directory as the sweep file
+- Create output in a `cabpack` folder one directory level up from the current directory
+- Automatically enable cabpack mode
+- Process all WAV files in the current directory (excluding the sweep file)
+- Use `selection.txt` and `ambient.txt` from the current directory if present
+
+**Perfect for:** Simply placing the tool in a folder with recorded WAV files and running it to create a complete cabpack in the parent directory's `cabpack` folder.
+
 ### Basic IR Generation
 
 Generate a standard impulse response from sweep and recorded files (any WAV format):
@@ -38,6 +60,131 @@ Generate a standard impulse response from sweep and recorded files (any WAV form
 ./valhallir-deconvolver --sweep sweep.wav --recorded recorded.wav --output ir.wav
 ```
 
+### Batch Processing (Directory Mode)
+
+Process all WAV files in a directory of recorded files. The tool will automatically detect if `--recorded` is a directory and process all WAV files found in it (including subdirectories):
+
+```sh
+./valhallir-deconvolver --sweep sweep.wav --recorded ./recordings/ --output ./ir_output/
+```
+
+This will:
+- Find all `.wav` files in the `./recordings/` directory (recursively)
+- Process each recorded file with the same sweep file
+- Generate IR files in the `./ir_output/` directory
+- Use the original recorded filename for each output IR (e.g., `recorded1.wav` → `recorded1.wav`)
+
+**Note:** If `--output` is a directory, it will be created automatically if it doesn't exist. If `--output` is a file path, single-file mode is used.
+
+**Example with options:**
+```sh
+./valhallir-deconvolver \
+  --sweep sweep.wav \
+  --recorded ./recordings/ \
+  --output ./ir_output/ \
+  --mpt \
+  --length-ms 50 \
+  --sample-rate 48000 \
+  --bit-depth 24
+```
+
+This processes all WAV files in `./recordings/` and generates both regular and MPT IRs in `./ir_output/`.
+
+### Cabpack Generation
+
+Generate a complete cabpack with IRs in multiple formats organized in a structured directory tree:
+
+**Simple mode (default - no options needed):**
+```sh
+./valhallir-deconvolver
+```
+This uses:
+- Current directory as recorded directory
+- `sweep.wav` in current directory as sweep file
+- `cabpack` folder one directory level up from current directory as output
+- Automatically enables cabpack mode
+- Excludes the sweep file from processing
+
+**Explicit mode:**
+```sh
+./valhallir-deconvolver \
+  --sweep sweep.wav \
+  --recorded ./recordings/ \
+  --output ./cabpack_output/ \
+  --cabpack
+```
+
+Both modes create a cabpack structure with:
+- **Format folders** named like `V2-1960STV 96000Hz-24bit 500ms` (extracted from WAV filenames)
+- **9 different formats** covering various sample rates, bit depths, and lengths:
+  - 44100Hz, 16bit, 170ms
+  - 44100Hz, 24bit, 170ms
+  - 44100Hz, 24bit, 500ms
+  - 48000Hz, 16bit, 170ms
+  - 48000Hz, 24bit, 170ms
+  - 48000Hz, 24bit, 500ms
+  - 48000Hz, 24bit, 1370ms
+  - 96000Hz, 24bit, 500ms
+  - 96000Hz, 24bit, 1370ms
+- **Subfolders** in each format: `close mics` (always), plus `ambient mics`, `selection`, and `mixes` (only if corresponding files/folders exist)
+- **RAW and MPT IRs** in `close mics/RAW` and `close mics/MPT` respectively
+- **IR visualizations** in the `plots` folder (96000Hz format)
+- **Automatic file organization** using `selection.txt` and `ambient.txt` files
+- **Mixes folder support**: Convert ready-to-use IRs from `mixes` folder to all formats
+
+The cabpack base name is automatically extracted from WAV filenames. For example, `V2-1960STV-d-SM7B-A1.wav` will create a cabpack named `V2-1960STV`.
+
+#### Automatic File Organization
+
+You can automatically populate the `selection` and `ambient mics` folders by placing text files in the recorded directory:
+
+- **`selection.txt`**: Lists filenames (one per line) to **copy** from `close mics` to `selection` folders. Both RAW and MPT versions are copied.
+- **`ambient.txt`**: Lists filenames (one per line) to **move** from `close mics` to `ambient mics` folders. Both RAW and MPT versions are moved.
+
+**Example `selection.txt`:**
+```
+V2-1960STV-d-SM7B-A1.wav
+V2-1960STV-d-SM7B-A2.wav
+V2-1960STV-d-SM7B-A3.wav
+```
+
+**Example `ambient.txt`:**
+```
+V2-1960STV-d-SM7B-B1.wav
+V2-1960STV-d-SM7B-B2.wav
+```
+
+**Notes:**
+- If `selection.txt` is missing, the `selection` folder is not created in any format folder.
+- If `ambient.txt` is missing, the `ambient mics` folder is not created in any format folder.
+- If the `mixes` folder is missing, the `mixes` subfolder is not created in any format folder.
+- Files are processed for all format folders automatically.
+- The `.wav` extension is automatically added if missing in the list files.
+
+#### Mixes Folder Support
+
+If a `mixes` folder exists in the recorded directory, it will be processed automatically:
+
+- **Ready-to-use IRs**: Place pre-processed IR files (already deconvolved) in the `mixes` folder
+- **Automatic conversion**: Each IR file is converted to all 9 cabpack formats
+- **Format storage**: Converted IRs are stored in the `mixes` subfolder of each format folder (only created if `mixes` folder exists)
+- **Plot generation**: IR visualizations are generated in the `plots/mixes` folder (only created if `mixes` folder exists)
+
+**Note:** If the `mixes` folder is missing in the recorded directory, the `mixes` subfolder is not created in any format folder, and no mix processing occurs.
+
+**Example structure:**
+```
+recorded/
+  ├── sweep.wav
+  ├── recorded1.wav
+  ├── recorded2.wav
+  └── mixes/
+      ├── mix1.wav
+      └── mix2.wav
+```
+
+This will create converted versions of `mix1.wav` and `mix2.wav` in all format folders under the `mixes` subfolder, plus plots in `plots/mixes/`.
+
 ### With Minimum Phase Transform
 
 Generate both regular and minimum phase IRs:
@@ -88,6 +235,33 @@ You can control the filter steepness (slope) with `--cut-slope` (in dB/octave, d
 
 This applies a 40 Hz low-cut and 18 kHz high-cut, both with a 24 dB/octave slope (steeper than the default 12).
 
+**Filter Implementation:** The tool uses Linkwitz-Riley design principles with Direct Form II Transposed filter structures for optimal numerical stability and minimal artifacts, even at very steep slopes (48+ dB/octave). The filters are designed to provide clean cutoffs without introducing unwanted artifacts in the frequency response.
+
+### Automatic Phase Correction
+
+Valhallir Deconvolver automatically detects and corrects phase-inverted recorded sweeps to ensure consistent IR polarity. This is especially important when mixing multiple IRs later.
+
+By default, the tool:
+- **Detects phase inversion** by analyzing the correlation between sweep and recorded signals
+- **Automatically corrects** phase-inverted recorded sweeps
+- **Ensures consistent polarity** across all generated IRs
+
+If you know your recorded sweep has the correct phase, you can disable automatic correction:
+
+```sh
+./valhallir-deconvolver --sweep sweep.wav --recorded recorded.wav --output ir.wav --no-phase-correction
+```
+
+### Forcing Phase Inversion (Testing/Manual Override)
+
+You can force the recorded sweep to be inverted (regardless of automatic detection) using:
+
+```sh
+./valhallir-deconvolver --sweep sweep.wav --recorded recorded.wav --output ir.wav --force-invert-phase
+```
+
+This is useful for testing or if you know your recorded sweep is out of phase and want to override the automatic detection.
+
 ### IR Visualization
 
 Generate frequency response and waveform plots of your IRs:
@@ -103,9 +277,11 @@ This creates:
 The plot includes:
 - **Frequency Response:** dB vs Hz with log frequency scale (20Hz-20kHz)
 - **Waveform:** Time-domain view of the first 10ms of the IR
-- **Valhallir Branding:** Logo and filename information
+- **Valhallir Branding:** Logo displayed in upper-left corner and filename information
 - **Professional Layout:** Clean, publication-ready visualization
 
+**Note:** The logo is embedded directly in the binary, so no external assets folder is required. The logo will always be available regardless of where the binary is located or how it's executed.
+
 ### Different Output Formats
 
 Generate IRs in different sample rates and bit depths:
@@ -156,9 +332,9 @@ Generate IRs in different sample rates and bit depths:
 
 | Flag | Description | Default | Required |
 |------|-------------|---------|----------|
-| `--sweep` | Path to sweep WAV file (any format) | - | Yes |
-| `--recorded` | Path to recorded WAV file (any format) | - | Yes |
-| `--output` | Path to output IR WAV file | - | Yes |
+| `--sweep` | Path to sweep WAV file (any format) | `sweep.wav` in current directory | No |
+| `--recorded` | Path to recorded WAV file or directory containing WAV files | Current directory | No |
+| `--output` | Path to output IR WAV file or directory for batch processing | `cabpack` folder one directory level up from recorded directory | No |
 | `--mpt` | Generate minimum phase transform IR | false | No |
 | `--sample-rate` | Output sample rate (44, 48, 88, 96 kHz) | 96000 | No |
 | `--bit-depth` | Output bit depth (16, 24, 32 bit) | 24 | No |
@@ -170,6 +346,9 @@ Generate IRs in different sample rates and bit depths:
 | `--highcut` | High-cut filter (low-pass) cutoff frequency in Hz (recorded sweep) | - | No |
 | `--cut-slope` | Filter slope in dB/octave (12, 24, 36, ...; default 12) | 12 | No |
 | `--plot-ir` | Generate frequency response and waveform plot | false | No |
+| `--no-phase-correction` | Disable automatic phase correction | false | No |
+| `--force-invert-phase` | Force inversion of the recorded sweep (manual override) | false | No |
+| `--cabpack` | Generate a cabpack with IRs in multiple formats organized in a directory tree | false | No |
 
 ## File Requirements
 
@@ -202,9 +381,11 @@ Generate IRs in different sample rates and bit depths:
 - You can change the fade duration with `--fade-ms`
 
 ### Filtering
-- You can apply a Butterworth low-cut (high-pass) and/or high-cut (low-pass) filter to the recorded sweep before IR extraction
+- You can apply a low-cut (high-pass) and/or high-cut (low-pass) filter to the recorded sweep before IR extraction
+- Uses Linkwitz-Riley design principles with Direct Form II Transposed structures for optimal stability
 - Use `--lowcut` and/or `--highcut` to specify cutoff frequencies in Hz
 - Use `--cut-slope` to control the filter steepness (12 dB/octave = gentle, 24+ = steeper)
+- The filter implementation is designed to provide clean cutoffs without artifacts, even at very steep slopes (48+ dB/octave)
 
 ### Deconvolution Process
 1. **FFT-based deconvolution** of recorded signal by sweep signal
@@ -225,6 +406,12 @@ Generate IRs in different sample rates and bit depths:
 - **Automatic File Naming:** Plots are saved with the same base name as the IR file
 - **High-Quality Output:** PNG format suitable for documentation and sharing
 
+### Phase Correction
+- **Automatic Detection:** Analyzes correlation between sweep and recorded signals to detect phase inversion
+- **Smart Correction:** Uses the first 100ms of signals for reliable phase analysis
+- **Consistent Polarity:** Ensures all IRs have the same phase polarity for easy mixing
+- **Optional Disable:** Use `--no-phase-correction` if you know the phase is correct
+
 ### Output Format Options
 - **Sample Rates:** 44.1kHz (CD), 48kHz (studio), 88.2kHz, 96kHz (high-res)
 - **Bit Depths:** 16-bit (CD), 24-bit (studio), 32-bit (high-res)
@@ -276,6 +463,18 @@ Generate IRs in different sample rates and bit depths:
   --mpt
 ```
 
+### Cabpack Generation
+```sh
+# Generate a complete cabpack with all formats
+./valhallir-deconvolver \
+  --sweep sweep.wav \
+  --recorded ./recordings/ \
+  --output ./cabpack_output/ \
+  --cabpack
+```
+
+This will process all WAV files in `./recordings/` and create a complete cabpack structure with IRs in 9 different formats, organized in folders with RAW and MPT versions, plus visualization plots.
+
 ## CI/CD & Multi-Platform Builds
 
 This project includes a Dagger pipeline for building binaries for multiple platforms:
@@ -354,4 +553,4 @@ The pipeline is defined in [`ci/dagger.go`](./ci/dagger.go). It outputs binaries
 
 ## Changelog
 
-See [CHANGELOG.md](./CHANGELOG.md) for version history and details. Current version: v1.0.0 
+See [CHANGELOG.md](./CHANGELOG.md) for version history and details. Current version: v1.2.0 

pkg/convolve/convolve.go

@@ -285,11 +285,25 @@ func realSlice(in []complex128) []float64 {
 }
 
 // Resample resamples audio data from one sample rate to another using linear interpolation
+// with proper anti-aliasing for downsampling
 func Resample(data []float64, fromSampleRate, toSampleRate int) []float64 {
 	if fromSampleRate == toSampleRate {
 		return data
 	}
 
+	// For downsampling, apply anti-aliasing filter to prevent aliasing
+	// Filter out frequencies above the target Nyquist frequency
+	if toSampleRate < fromSampleRate {
+		// Calculate target Nyquist frequency (slightly below to ensure clean cutoff)
+		targetNyquist := float64(toSampleRate) / 2.0
+		// Use 90% of Nyquist as cutoff to ensure clean anti-aliasing
+		antiAliasCutoff := targetNyquist * 0.9
+
+		// Apply steep low-pass filter to prevent aliasing
+		// Use 48 dB/octave slope for clean anti-aliasing
+		data = CascadeHighcut(data, fromSampleRate, antiAliasCutoff, 48)
+	}
+
 	// Calculate the resampling ratio
 	ratio := float64(toSampleRate) / float64(fromSampleRate)
 	newLength := int(float64(len(data)) * ratio)
@@ -342,13 +356,71 @@ func FadeOutLinear(data []float64, fadeSamples int) []float64 {
 	return out
 }
 
+// biquadFilterState holds the state for a biquad filter
+type biquadFilterState struct {
+	x1, x2 float64 // input history (intermediate values in DF2T)
+}
+
+// applyBiquadDF2T applies a biquad filter using Direct Form II Transposed (more numerically stable)
+// This form is preferred for cascaded filters as it reduces numerical errors
+// Uses high-precision calculations and proper state management for maximum quality
+func applyBiquadDF2T(data []float64, b0, b1, b2, a1, a2 float64, state *biquadFilterState) []float64 {
+	out := make([]float64, len(data))
+
+	// Initialize state if nil
+	if state == nil {
+		state = &biquadFilterState{}
+	}
+
+	// Direct Form II Transposed implementation with high precision
+	// This structure minimizes numerical errors in cascaded filters
+	for i := 0; i < len(data); i++ {
+		x := data[i]
+		// Compute intermediate value (w[n] = x[n] - a1*w[n-1] - a2*w[n-2])
+		w := x - a1*state.x1 - a2*state.x2
+		// Compute output (y[n] = b0*w[n] + b1*w[n-1] + b2*w[n-2])
+		y := b0*w + b1*state.x1 + b2*state.x2
+		out[i] = y
+		// Update state (shift delay line)
+		state.x2 = state.x1
+		state.x1 = w
+	}
+
+	return out
+}
+
+// warmupFilter applies a filter with warm-up to avoid transients
+// This is especially important for high-quality filtering
+func warmupFilter(data []float64, b0, b1, b2, a1, a2 float64) []float64 {
+	if len(data) == 0 {
+		return data
+	}
+
+	// Use first sample value for warm-up to avoid transients
+	warmupValue := data[0]
+	state := &biquadFilterState{}
+
+	// Warm up the filter with a few samples of the first value
+	// This initializes the filter state properly
+	warmupSamples := 10
+	warmupData := make([]float64, warmupSamples)
+	for i := range warmupData {
+		warmupData[i] = warmupValue
+	}
+	_ = applyBiquadDF2T(warmupData, b0, b1, b2, a1, a2, state)
+
+	// Now apply filter to actual data with warmed-up state
+	return applyBiquadDF2T(data, b0, b1, b2, a1, a2, state)
+}
+
 // ApplyLowpassButterworth applies a 2nd-order Butterworth low-pass filter to the data.
 // cutoffHz: cutoff frequency in Hz, sampleRate: sample rate in Hz.
+// Uses Direct Form II Transposed for better numerical stability.
 func ApplyLowpassButterworth(data []float64, sampleRate int, cutoffHz float64) []float64 {
 	if cutoffHz <= 0 || cutoffHz >= float64(sampleRate)/2 {
 		return data
 	}
-	// Biquad coefficients
+	// Biquad coefficients for Butterworth low-pass
 	w0 := 2 * math.Pi * cutoffHz / float64(sampleRate)
 	cosw0 := math.Cos(w0)
 	sinw0 := math.Sin(w0)
@@ -362,35 +434,26 @@ func ApplyLowpassButterworth(data []float64, sampleRate int, cutoffHz float64) [
 	a1 := -2 * cosw0
 	a2 := 1 - alpha
 
-	// Normalize
+	// Normalize coefficients
 	b0 /= a0
 	b1 /= a0
 	b2 /= a0
 	a1 /= a0
 	a2 /= a0
 
-	// Apply filter (Direct Form I)
-	out := make([]float64, len(data))
-	var x1, x2, y1, y2 float64
-	for i := 0; i < len(data); i++ {
-		x0 := data[i]
-		y0 := b0*x0 + b1*x1 + b2*x2 - a1*y1 - a2*y2
-		out[i] = y0
-		x2 = x1
-		x1 = x0
-		y2 = y1
-		y1 = y0
-	}
-	return out
+	// Apply filter using Direct Form II Transposed
+	state := &biquadFilterState{}
+	return applyBiquadDF2T(data, b0, b1, b2, a1, a2, state)
 }
 
 // ApplyHighpassButterworth applies a 2nd-order Butterworth high-pass filter to the data.
 // cutoffHz: cutoff frequency in Hz, sampleRate: sample rate in Hz.
+// Uses Direct Form II Transposed for better numerical stability.
 func ApplyHighpassButterworth(data []float64, sampleRate int, cutoffHz float64) []float64 {
 	if cutoffHz <= 0 || cutoffHz >= float64(sampleRate)/2 {
 		return data
 	}
-	// Biquad coefficients
+	// Biquad coefficients for Butterworth high-pass
 	w0 := 2 * math.Pi * cutoffHz / float64(sampleRate)
 	cosw0 := math.Cos(w0)
 	sinw0 := math.Sin(w0)
@@ -404,52 +467,235 @@ func ApplyHighpassButterworth(data []float64, sampleRate int, cutoffHz float64)
 	a1 := -2 * cosw0
 	a2 := 1 - alpha
 
-	// Normalize
+	// Normalize coefficients
 	b0 /= a0
 	b1 /= a0
 	b2 /= a0
 	a1 /= a0
 	a2 /= a0
 
-	// Apply filter (Direct Form I)
-	out := make([]float64, len(data))
-	var x1, x2, y1, y2 float64
-	for i := 0; i < len(data); i++ {
-		x0 := data[i]
-		y0 := b0*x0 + b1*x1 + b2*x2 - a1*y1 - a2*y2
-		out[i] = y0
-		x2 = x1
-		x1 = x0
-		y2 = y1
-		y1 = y0
-	}
-	return out
+	// Apply filter using Direct Form II Transposed
+	state := &biquadFilterState{}
+	return applyBiquadDF2T(data, b0, b1, b2, a1, a2, state)
 }
 
 // CascadeLowcut applies the low-cut (high-pass) filter multiple times for steeper slopes.
-// slopeDb: 12, 24, 36, ... (dB/octave)
+// Uses Linkwitz-Riley design principles with high-quality implementation for maximum precision.
+// Features: proper warm-up, high-precision coefficients, optimized cascade structure.
+// slopeDb: 12, 24, 36, 48, ... (dB/octave)
 func CascadeLowcut(data []float64, sampleRate int, cutoffHz float64, slopeDb int) []float64 {
 	if slopeDb < 12 {
 		slopeDb = 12
 	}
-	n := slopeDb / 12
-	out := data
-	for i := 0; i < n; i++ {
-		out = ApplyHighpassButterworth(out, sampleRate, cutoffHz)
+	if cutoffHz <= 0 || cutoffHz >= float64(sampleRate)/2 {
+		return data
 	}
+
+	n := slopeDb / 12
+	if n == 0 {
+		return data
+	}
+
+	// High-quality filter implementation with proper coefficient calculation
+	// Use high precision for coefficient calculation
+	w0 := 2.0 * math.Pi * cutoffHz / float64(sampleRate)
+
+	// Pre-calculate trigonometric functions once for precision
+	cosw0 := math.Cos(w0)
+	sinw0 := math.Sin(w0)
+
+	// Butterworth Q factor (1/sqrt(2) ≈ 0.7071067811865476 for maximally flat response)
+	const butterworthQ = 0.7071067811865476
+	alpha := sinw0 / (2.0 * butterworthQ)
+
+	// High-pass Butterworth coefficients
+	b0 := (1.0 + cosw0) / 2.0
+	b1 := -(1.0 + cosw0)
+	b2 := (1.0 + cosw0) / 2.0
+	a0 := 1.0 + alpha
+	a1 := -2.0 * cosw0
+	a2 := 1.0 - alpha
+
+	// Normalize coefficients for Direct Form II Transposed
+	// This ensures proper scaling and numerical stability
+	b0 /= a0
+	b1 /= a0
+	b2 /= a0
+	a1 /= a0
+	a2 /= a0
+
+	// Apply cascaded filters with proper warm-up for each stage
+	// This ensures clean filtering without transients, especially important for MPT
+	out := make([]float64, len(data))
+	copy(out, data)
+
+	for stage := 0; stage < n; stage++ {
+		// Use warm-up for first stage to avoid transients
+		// Subsequent stages benefit from the already-filtered signal
+		if stage == 0 {
+			out = warmupFilter(out, b0, b1, b2, a1, a2)
+		} else {
+			// For subsequent stages, use fresh state but no warm-up needed
+			// as the signal is already filtered
+			state := &biquadFilterState{}
+			filtered := applyBiquadDF2T(out, b0, b1, b2, a1, a2, state)
+			copy(out, filtered)
+		}
+	}
+
 	return out
 }
 
 // CascadeHighcut applies the high-cut (low-pass) filter multiple times for steeper slopes.
-// slopeDb: 12, 24, 36, ... (dB/octave)
+// Uses Linkwitz-Riley design principles with high-quality implementation for maximum precision.
+// Features: proper warm-up, high-precision coefficients, optimized cascade structure.
+// slopeDb: 12, 24, 36, 48, ... (dB/octave)
 func CascadeHighcut(data []float64, sampleRate int, cutoffHz float64, slopeDb int) []float64 {
 	if slopeDb < 12 {
 		slopeDb = 12
 	}
-	n := slopeDb / 12
-	out := data
-	for i := 0; i < n; i++ {
-		out = ApplyLowpassButterworth(out, sampleRate, cutoffHz)
+	if cutoffHz <= 0 || cutoffHz >= float64(sampleRate)/2 {
+		return data
 	}
+
+	n := slopeDb / 12
+	if n == 0 {
+		return data
+	}
+
+	// High-quality filter implementation with proper coefficient calculation
+	// Use high precision for coefficient calculation
+	w0 := 2.0 * math.Pi * cutoffHz / float64(sampleRate)
+
+	// Pre-calculate trigonometric functions once for precision
+	cosw0 := math.Cos(w0)
+	sinw0 := math.Sin(w0)
+
+	// Butterworth Q factor (1/sqrt(2) ≈ 0.7071067811865476 for maximally flat response)
+	const butterworthQ = 0.7071067811865476
+	alpha := sinw0 / (2.0 * butterworthQ)
+
+	// Low-pass Butterworth coefficients
+	b0 := (1.0 - cosw0) / 2.0
+	b1 := 1.0 - cosw0
+	b2 := (1.0 - cosw0) / 2.0
+	a0 := 1.0 + alpha
+	a1 := -2.0 * cosw0
+	a2 := 1.0 - alpha
+
+	// Normalize coefficients for Direct Form II Transposed
+	// This ensures proper scaling and numerical stability
+	b0 /= a0
+	b1 /= a0
+	b2 /= a0
+	a1 /= a0
+	a2 /= a0
+
+	// Apply cascaded filters with proper warm-up for each stage
+	// This ensures clean filtering without transients, especially important for MPT
+	out := make([]float64, len(data))
+	copy(out, data)
+
+	for stage := 0; stage < n; stage++ {
+		// Use warm-up for first stage to avoid transients
+		// Subsequent stages benefit from the already-filtered signal
+		if stage == 0 {
+			out = warmupFilter(out, b0, b1, b2, a1, a2)
+		} else {
+			// For subsequent stages, use fresh state but no warm-up needed
+			// as the signal is already filtered
+			state := &biquadFilterState{}
+			filtered := applyBiquadDF2T(out, b0, b1, b2, a1, a2, state)
+			copy(out, filtered)
+		}
+	}
+
 	return out
 }
+
+// min returns the minimum of three integers
+func min(a, b, c int) int {
+	if a <= b && a <= c {
+		return a
+	}
+	if b <= a && b <= c {
+		return b
+	}
+	return c
+}
+
+// DetectPhaseInversion detects if the recorded sweep is phase-inverted compared to the sweep
+// by computing the normalized cross-correlation over a range of lags
+func DetectPhaseInversion(sweep, recorded []float64) bool {
+	if len(sweep) == 0 || len(recorded) == 0 {
+		return false
+	}
+
+	windowSize := min(len(sweep), len(recorded), 9600) // 100ms at 96kHz
+	sweepWindow := sweep[:windowSize]
+	recordedWindow := recorded[:windowSize]
+
+	maxLag := 500 // +/- 500 samples (~5ms)
+	bestCorr := 0.0
+	bestLag := 0
+
+	for lag := -maxLag; lag <= maxLag; lag++ {
+		var corr, sweepSum, recordedSum, sweepSumSq, recordedSumSq float64
+		count := 0
+		for i := 0; i < windowSize; i++ {
+			j := i + lag
+			if j < 0 || j >= windowSize {
+				continue
+			}
+			corr += sweepWindow[i] * recordedWindow[j]
+			sweepSum += sweepWindow[i]
+			recordedSum += recordedWindow[j]
+			sweepSumSq += sweepWindow[i] * sweepWindow[i]
+			recordedSumSq += recordedWindow[j] * recordedWindow[j]
+			count++
+		}
+		if count == 0 {
+			continue
+		}
+		sweepMean := sweepSum / float64(count)
+		recordedMean := recordedSum / float64(count)
+		sweepVar := sweepSumSq/float64(count) - sweepMean*sweepMean
+		recordedVar := recordedSumSq/float64(count) - recordedMean*recordedMean
+		if sweepVar <= 0 || recordedVar <= 0 {
+			continue
+		}
+		corrCoeff := (corr/float64(count) - sweepMean*recordedMean) / math.Sqrt(sweepVar*recordedVar)
+		if math.Abs(corrCoeff) > math.Abs(bestCorr) {
+			bestCorr = corrCoeff
+			bestLag = lag
+		}
+	}
+
+	log.Printf("[deconvolve] Phase cross-correlation: best lag = %d, coeff = %.4f", bestLag, bestCorr)
+	return bestCorr < 0.0
+}
+
+// InvertPhase inverts the phase of the audio data by negating all samples
+func InvertPhase(data []float64) []float64 {
+	inverted := make([]float64, len(data))
+	for i, sample := range data {
+		inverted[i] = -sample
+	}
+	return inverted
+}
+
+// DeconvolveWithPhaseCorrection extracts the impulse response with automatic phase correction
+func DeconvolveWithPhaseCorrection(sweep, recorded []float64) []float64 {
+	// Detect if recorded sweep is phase-inverted
+	isInverted := DetectPhaseInversion(sweep, recorded)
+
+	if isInverted {
+		log.Printf("[deconvolve] Detected phase inversion in recorded sweep, correcting...")
+		recorded = InvertPhase(recorded)
+	} else {
+		log.Printf("[deconvolve] Phase alignment verified")
+	}
+
+	// Perform normal deconvolution
+	return Deconvolve(sweep, recorded)
+}

pkg/plot/plot.go

@@ -1,17 +1,16 @@
 package plot
 
 import (
+	"bytes"
 	"fmt"
+	"image/color"
+	"image/png"
 	"math"
 	"math/cmplx"
 	"os"
 	"path/filepath"
 	"strings"
 
-	"image/png"
-
-	"image/color"
-
 	"github.com/mjibson/go-dsp/fft"
 	"gonum.org/v1/plot"
 	"gonum.org/v1/plot/font"
@@ -22,7 +21,8 @@ import (
 )
 
 // PlotIR plots the frequency response (magnitude in dB vs. frequency in Hz) of the IR to a PNG file
-func PlotIR(ir []float64, sampleRate int, irFileName string) error {
+// logoData is the embedded logo image data (can be nil if not available)
+func PlotIR(ir []float64, sampleRate int, irFileName string, logoData []byte) error {
 	if len(ir) == 0 {
 		return nil
 	}
@@ -63,7 +63,7 @@ func PlotIR(ir []float64, sampleRate int, irFileName string) error {
 	}
 	fmt.Printf("[PlotIR] minDb in plotted range: %.2f dB at %.2f Hz\n", minDb, minDbFreq)
 	p := plot.New()
-	p.Title.Text = "IR Frequency Response (dB, 2048-sample window)"
+	p.Title.Text = "IR Frequency Response"
 	p.X.Label.Text = "Frequency (Hz)"
 	p.Y.Label.Text = "Magnitude (dB)"
 	p.X.Scale = plot.LogScale{}
@@ -111,7 +111,7 @@ func PlotIR(ir []float64, sampleRate int, irFileName string) error {
 
 	// --- Time-aligned waveform plot ---
 	p2 := plot.New()
-	p2.Title.Text = "IR Waveform (Time Aligned)"
+	p2.Title.Text = "IR Waveform"
 	p2.X.Label.Text = "Time (ms)"
 	p2.Y.Label.Text = "Amplitude"
 	// Prepare waveform data (only first 10ms)
@@ -141,16 +141,14 @@ func PlotIR(ir []float64, sampleRate int, irFileName string) error {
 	dc := draw.New(img)
 
 	// Draw logo at the top left, headline to the right, IR filename below
-	logoPath := "assets/logo.png"
+	// Logo is embedded in the binary
 	logoW := 2.4 * vg.Inch  // doubled size
 	logoH := 0.68 * vg.Inch // doubled size
 	logoX := 0.3 * vg.Inch
 	logoY := height + 0.2*vg.Inch // move logo down by an additional ~10px
 	logoDrawn := false
-	f, err := os.Open(logoPath)
-	if err == nil {
-		defer f.Close()
-		logoImg, err := png.Decode(f)
+	if len(logoData) > 0 {
+		logoImg, err := png.Decode(bytes.NewReader(logoData))
 		if err == nil {
 			rect := vg.Rectangle{
 				Min: vg.Point{X: logoX, Y: logoY},
@@ -201,11 +199,11 @@ func PlotIR(ir []float64, sampleRate int, irFileName string) error {
 	irNameWithoutExt := strings.TrimSuffix(irBase, filepath.Ext(irBase))
 	plotFileName := filepath.Join(irDir, irNameWithoutExt+".png")
 
-	f, err = os.Create(plotFileName)
+	plotFile, err := os.Create(plotFileName)
 	if err != nil {
 		return err
 	}
-	defer f.Close()
-	_, err = vgimg.PngCanvas{Canvas: img}.WriteTo(f)
+	defer plotFile.Close()
+	_, err = vgimg.PngCanvas{Canvas: img}.WriteTo(plotFile)
 	return err
 }

pkg/wav/reader.go

@@ -36,6 +36,15 @@ func toMono(data []float64, channels int) []float64 {
 
 // ReadWAVFile reads a WAV file and returns its PCM data as float64 (resampled to 96kHz mono)
 func ReadWAVFile(filePath string) (*WAVData, error) {
+	// Check if path is a directory
+	info, err := os.Stat(filePath)
+	if err != nil {
+		return nil, fmt.Errorf("failed to access file %s: %w", filePath, err)
+	}
+	if info.IsDir() {
+		return nil, fmt.Errorf("path %s is a directory, not a WAV file", filePath)
+	}
+
 	file, err := os.Open(filePath)
 	if err != nil {
 		return nil, fmt.Errorf("failed to open file %s: %w", filePath, err)