File src/remiaudio/dsp/reverbs/zita-paraeq.cr from the latest check-in


     1
     2
     3
     4
     5
     6
     7
     8
     9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    56
    57
    58
    59
    60
    61
    62
    63
    64
    65
    66
    67
    68
    69
    70
    71
    72
    73
    74
    75
    76
    77
    78
    79
    80
    81
    82
    83
    84
    85
    86
    87
    88
    89
    90
    91
    92
    93
    94
    95
    96
    97
    98
    99
   100
   101
   102
   103
   104
   105
   106
   107
   108
   109
   110
   111
   112
   113
   114
   115
   116
   117
   118
   119
   120
   121
   122
   123
   124
   125
   126
   127
   128
   129
   130
   131
   132
   133
   134
   135
   136
   137
   138
   139
   140
   141
   142
   143
   144
   145
   146
   147
   148
   149
   150
   151
   152
   153
   154
   155
   156
   157
   158
   159
   160
   161
   162
   163
   164
   165
   166
   167
   168
   169
   170
   171
   172
   173
   174
   175
   176
   177
   178
   179
   180
   181
   182
   183
   184
   185
   186
   187
   188
   189
   190
   191
   192
   193
   194
   195
   196
   197
   198
   199
   200
   201
   202
   203
   204
   205
   206
   207
   208
   209
   210
   211
   212
   213
   214
   215
   216
   217
   218
   219
   220
   221
   222
   223
   224
   225
   226
   227
   228
   229
   230
   231
   232
   233
   234
   235
   236
   237
   238
   239
   240
   241
   242
   243
   244
   245
   246
   247
   248
   249
   250
   251
   252
   253
   254
   255
   256
   257
   258
   259
   260
   261
   262
   263
   264
   265
   266
   267
   268
   269
   270
   271
   272
   273
   274
   275
   276
   277
   278
   279
   280
   281
   282
   283
   284
   285
   286
   287
   288
   289
   290
   291
   292
   293
   294
   295
   296
   297
   298
   299
   300
   301
   302
   303
   304
   305
   306
   307
   308
   309
   310
   311
   312
   313
   314
   315
   316
   317
   318
   319
   320
   321
   322
   323
   324
   325
   326
   327
   328
   329
   330
   331
   332
   333
   334
   335
   336
   337
   338
   339
   340
   341
   342
   343
   344
   345
   346
   347
   348
   349
   350
   351
   352
   353
   354
   355
   356
   357
   358
   359
   360
   361
   362
   363
   364
   365
   366
   367
   368
   369
   370
   371
   372
   373
   374
   375
   376
   377
   378
   379
   380
   381
   382
   383
   384
   385
   386
   387
   388
   389
   390
   391
   392
   393
   394
   395
   396
   397
   398
   399
   400
   401
   402
   403
   404
   405
   406
   407
   408
   409
   410
   411
   412
   413
   414
   415
   416
   417
   418
   419
   420
   421
   422
   423
   424
   425
   426
   427
   428
   429
   430
   431
   432
   433
   434
   435
   436
   437
   438
   439
   440
   441
   442
   443
   444
   445
   446
   447
   448
   449
   450
   451
   452
   453
   454
   455
   456
   457
   458
   459
   460
   461
   462
   463
   464
   465
   466
   467
   468
   469
   470
   471
   472
   473
   474
   475
   476
   477
   478
   479
   480
   481
   482
   483
#### RemiAudio
#### Copyright (C) 2022-2024 Remilia Scarlet <remilia@posteo.jp>
####   Based on Zita-Rev1
####   Copyright (C) 2003-2017 Fons Adriaensen <fons@linuxaudio.org>
####
#### This program is free software; you can redistribute it and/or modify it
#### under the terms of the GNU General Public License as published by the Free
#### Software Foundation; either version 3 of the License, or (at your option)
#### any later version.
####
#### This program is distributed in the hope that it will be useful, but WITHOUT
#### ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
#### FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for
#### more details.
####
#### You should have received a copy of the GNU General Public License along
#### with this program.  If not, see <http://www.gnu.org/licenses/>.
require "./reverb"

module RemiAudio::DSP
  class ZitaReverb < Reverb
    # A parametric equalizer specifically for `ZitaReverb` and
    # `ZitaReverbSimple`.
    private class ParametricEQ
      enum Mode
        Bypass = 1
        Static
        Smooth
        MaxCh
      end

      @touch0 : Int16 = 0
      @touch1 : Int16 = 0
      @state  : Mode = Mode::Bypass
      @bypass : Bool = false
      @sampleRate : Float64 = 0.0
      @invSampleRate : Float64 = 0.0
      getter g0 : Float64 = 1.0
      getter g1 : Float64 = 1.0
      getter f0 : Float64 = 1.0e3
      getter f1 : Float64 = 1.0e3
      @c1 : Float64 = 0.0
      @dc1 : Float64 = 0.0
      @c2 : Float64 = 0.0
      @dc2 : Float64 = 0.0
      @gg : Float64 = 0.0
      @dgg : Float64 = 0.0
      @z1 : Array(Float64) = Array.new(PARA_EQ_MAX_CH, 0.0f64)
      @z2 : Array(Float64) = Array.new(PARA_EQ_MAX_CH, 0.0f64)

      def initialize
        self.sampleRate = 44100.0
      end

      @[AlwaysInline]
      def sampleRate=(@sampleRate : Float64) : Nil
        @invSampleRate = 1.0 / @sampleRate
        reset
      end

      @[AlwaysInline]
      def setParam(freq : Float64, gain : Float64) : Nil
        @f0 = freq
        @g0 = 10.0f64 ** (0.05 * gain)
        @touch0 += 1
      end

      @[AlwaysInline]
      def reset : Nil
        @z1.fill(0.0f64)
        @z2.fill(0.0f64)
      end

      def prepare(numFrames : Int) : Nil
        update : Bool = false
        g : Float64 = 0.0
        f : Float64 = 0.0

        unless @touch0 == @touch1
          g = @g0
          f = @f0

          unless g == @g1
            update = true
            if g > (@g1 + @g1)
              @g1 *= 2.0
            elsif @g1 > (g + g)
              @g1 *= 0.5
            else
              @g1 = g
            end
          end

          unless f == @f1
            update = true
            if f > (@f1 + @f1)
              @f1 *= 2.0
            elsif @f1 > (f + f)
              @f1 *= 0.5
            else
              @f1 = f
            end
          end

          if update
            if @state.bypass? && @g1 == 1.0
              calculate(0, @g1, @f1)
            else
              @state = Mode::Smooth
              calculate(numFrames, @g1, @f1)
            end
          else
            @touch1 = @touch0
            if (@g1 - 1).abs < 0.001
              @state = Mode::Bypass
              reset
            else
              @state = Mode::Static
            end
          end
        end
      end

      @[AlwaysInline]
      def calculate(numFrames : Int, gain : Float64, freq : Float64) : Nil
        freq = freq * Math::PI * @invSampleRate
        b : Float64 = (2 * freq) / Math.sqrt(gain)
        gg : Float64 = 0.5 * (gain - 1.0)
        c1 : Float64 = -RemiMath.fastCos(2 * freq)
        c2 : Float64 = (1 - b) / (1 + b)

        if numFrames > 0
          @dc1 = (c1 - @c1) / numFrames + 1.0e-30f64
          @dc2 = (c2 - @c2) / numFrames + 1.0e-30f64
          @dgg = (gg - @gg) / numFrames + 1.0e-30f64
        else
          @c1 = c1
          @c2 = c2
          @gg = gg
        end
      end

      private macro defineProcessFuncs(numType, fromDataFn, toDataFn)
        {% if (toDataFn && !fromDataFn) || (!toDataFn && fromDataFn) %}
          {% raise "fromDataFn and toDataFn must both be specified, or neither" %}
        {% end %}

        private def processSmooth(dataLeft : Array({{numType}})|Slice({{numType}}),
                                  dataRight : Array({{numType}})|Slice({{numType}})) : Nil
          c1 : Float64 = 0.0
          c2 : Float64 = 0.0
          gg : Float64 = 0.0
          z1 : Float64 = 0.0
          z2 : Float64 = 0.0
          x  : Float64 = 0.0
          y  : Float64 = 0.0

          # Left...
          z1 = @z1.unsafe_fetch(0)
          z2 = @z2.unsafe_fetch(0)
          c1 = @c1
          c2 = @c2
          gg = @gg

          dataLeft.size.times do |j|
            c1 += @dc1
            c2 += @dc2
            gg += @dgg

            {% if fromDataFn %}
              x = dataLeft.unsafe_fetch(j).{{fromDataFn.id}}
            {% else %}
              x = dataLeft.unsafe_fetch(j)
            {% end %}

            y = x - (c2 * z2)

            {% if toDataFn %}
              dataLeft.unsafe_put(j, (x - gg * (z2 + c2 * y - x)).{{toDataFn.id}})
            {% else %}
              dataLeft.unsafe_put(j, x - gg * (z2 + c2 * y - x))
            {% end %}

            y -= c1 * z1
            z2 = z1 + c1 * y
            z1 = y + 1.0e-20f64
          end
          @z1.unsafe_put(0, z1)
          @z2.unsafe_put(0, z2)

          # Right...
          z1 = @z1.unsafe_fetch(1)
          z2 = @z2.unsafe_fetch(1)
          c1 = @c1
          c2 = @c2
          gg = @gg
          dataRight.size.times do |j|
            c1 += @dc1
            c2 += @dc2
            gg += @dgg

            {% if fromDataFn %}
              x = dataRight.unsafe_fetch(j).{{fromDataFn.id}}
            {% else %}
              x = dataRight.unsafe_fetch(j)
            {% end %}

            x = dataRight[j]
            y = x - c2 * z2

            {% if toDataFn %}
              dataRight.unsafe_put(j, (x - gg * (z2 + c2 * y - x)).{{toDataFn.id}})
            {% else %}
              dataRight.unsafe_put(j, x - gg * (z2 + c2 * y - x))
            {% end %}

            y -= c1 * z1
            z2 = z1 + c1 * y
            z1 = y + 1.0e-20f64
          end
          @z1.unsafe_put(1, z1)
          @z2.unsafe_put(1, z2)

          # Finish up
          @c1 = c1
          @c2 = c2
          @gg = gg
        end

        private def processNonSmooth(dataLeft : Array({{numType}})|Slice({{numType}}),
                                     dataRight : Array({{numType}})|Slice({{numType}})) : Nil
          c1 : Float64 = 0.0
          c2 : Float64 = 0.0
          gg : Float64 = 0.0
          z1 : Float64 = 0.0
          z2 : Float64 = 0.0
          x  : Float64 = 0.0
          y  : Float64 = 0.0

          # Left...
          z1 = @z1.unsafe_fetch(0)
          z2 = @z2.unsafe_fetch(0)

          dataLeft.size.times do |j|
            {% if fromDataFn %}
              x = dataLeft.unsafe_fetch(j).{{fromDataFn.id}}
            {% else %}
              x = dataLeft.unsafe_fetch(j)
            {% end %}

            y = x - (c2 * z2)

            {% if toDataFn %}
              dataLeft.unsafe_put(j, (x - gg * (z2 + c2 * y - x)).{{toDataFn.id}})
            {% else %}
              dataLeft.unsafe_put(j, x - gg * (z2 + c2 * y - x))
            {% end %}

            y -= c1 * z1
            z2 = z1 + (c1 * y)
            z1 = y + 1.0e-20f64
          end
          @z1.unsafe_put(0, z1)
          @z2.unsafe_put(0, z2)

          # Right...
          z1 = @z1.unsafe_fetch(1)
          z2 = @z2.unsafe_fetch(1)
          dataRight.size.times do |j|
            {% if fromDataFn %}
              x = dataRight.unsafe_fetch(j).{{fromDataFn.id}}
            {% else %}
              x = dataRight.unsafe_fetch(j)
            {% end %}

            y = x - (c2 * z2)
            {% if toDataFn %}
              dataRight.unsafe_put(j, (x - gg * (z2 + c2 * y - x)).{{toDataFn.id}})
            {% else %}
              dataRight.unsafe_put(j, x - gg * (z2 + c2 * y - x))
            {% end %}

            y -= (c1 * z1)
            z2 = z1 + (c1 * y)
            z1 = y + 1.0e-20f64
          end
          @z1.unsafe_put(1, z1)
          @z2.unsafe_put(1, z2)
        end

        # Interleaved version
        private def processSmooth(data : Array({{numType}})|Slice({{numType}})) : Nil
          c1 : Float64 = 0.0
          c2 : Float64 = 0.0
          gg : Float64 = 0.0
          z1 : Float64 = 0.0
          z2 : Float64 = 0.0
          x  : Float64 = 0.0
          y  : Float64 = 0.0
          j = 0

          # Left...
          z1 = @z1.unsafe_fetch(0)
          z2 = @z2.unsafe_fetch(0)
          c1 = @c1
          c2 = @c2
          gg = @gg
          while j < data.size
            c1 += @dc1
            c2 += @dc2
            gg += @dgg

            {% if fromDataFn %}
              x = data.unsafe_fetch(j).{{fromDataFn.id}}
            {% else %}
              x = data.unsafe_fetch(j)
            {% end %}

            y = x - (c2 * z2)

            {% if toDataFn %}
              data.unsafe_put(j, (x - gg * (z2 + c2 * y - x)).{{toDataFn.id}})
            {% else %}
              data.unsafe_put(j, x - gg * (z2 + c2 * y - x))
            {% end %}

            y -= c1 * z1
            z2 = z1 + c1 * y
            z1 = y + 1.0e-20f64
            j += 2
          end
          @z1.unsafe_put(0, z1)
          @z2.unsafe_put(0, z2)

          # Right...
          z1 = @z1.unsafe_fetch(1)
          z2 = @z2.unsafe_fetch(1)
          c1 = @c1
          c2 = @c2
          gg = @gg
          j = 1
          while j < data.size
            c1 += @dc1
            c2 += @dc2
            gg += @dgg

            {% if fromDataFn %}
              x = data.unsafe_fetch(j).{{fromDataFn.id}}
            {% else %}
              x = data.unsafe_fetch(j)
            {% end %}

            y = x - c2 * z2

            {% if toDataFn %}
              data.unsafe_put(j, (x - gg * (z2 + c2 * y - x)).{{toDataFn.id}})
            {% else %}
              data.unsafe_put(j, x - gg * (z2 + c2 * y - x))
            {% end %}

            y -= c1 * z1
            z2 = z1 + c1 * y
            z1 = y + 1.0e-20f64
            j += 2
          end
          @z1.unsafe_put(1, z1)
          @z2.unsafe_put(1, z2)

          # Finish up
          @c1 = c1
          @c2 = c2
          @gg = gg
        end

        # Interleaved version
        private def processNonSmooth(data : Array({{numType}})|Slice({{numType}})) : Nil
          c1 : Float64 = 0.0
          c2 : Float64 = 0.0
          gg : Float64 = 0.0
          z1 : Float64 = 0.0
          z2 : Float64 = 0.0
          x  : Float64 = 0.0
          y  : Float64 = 0.0
          j = 0

          # Left...
          z1 = @z1.unsafe_fetch(0)
          z2 = @z2.unsafe_fetch(0)
          while j < data.size
            {% if fromDataFn %}
              x = data.unsafe_fetch(j).{{fromDataFn.id}}
            {% else %}
              x = data.unsafe_fetch(j)
            {% end %}

            y = x - (c2 * z2)

            {% if toDataFn %}
              data.unsafe_put(j, (x - gg * (z2 + c2 * y - x)).{{toDataFn.id}})
            {% else %}
              data.unsafe_put(j, x - gg * (z2 + c2 * y - x))
            {% end %}

            y -= c1 * z1
            z2 = z1 + (c1 * y)
            z1 = y + 1.0e-20f64
            j += 2
          end
          @z1.unsafe_put(0, z1)
          @z2.unsafe_put(0, z2)

          # Right...
          z1 = @z1.unsafe_fetch(1)
          z2 = @z2.unsafe_fetch(1)
          j = 1
          while j < data.size
            {% if fromDataFn %}
              x = data.unsafe_fetch(j).{{fromDataFn.id}}
            {% else %}
              x = data.unsafe_fetch(j)
            {% end %}

            y = x - (c2 * z2)

            {% if toDataFn %}
              data.unsafe_put(j, (x - gg * (z2 + c2 * y - x)).{{toDataFn.id}})
            {% else %}
              data.unsafe_put(j, x - gg * (z2 + c2 * y - x))
            {% end %}

            y -= (c1 * z1)
            z2 = z1 + (c1 * y)
            z1 = y + 1.0e-20f64
            j += 2
          end
          @z1.unsafe_put(1, z1)
          @z2.unsafe_put(1, z2)
        end

      end

      defineProcessFuncs(Float64, nil, nil)
      defineProcessFuncs(Float32, :to_f64!, :to_f32!)

      {% begin %}
        {% for typ in [:Float64, :Float32] %}
          @[AlwaysInline]
          def process(dataLeft : Array({{typ.id}})|Slice({{typ.id}}),
                      dataRight : Array({{typ.id}})|Slice({{typ.id}})) : Nil
            if @state.smooth?
              processSmooth(dataLeft, dataRight)
            else
              processNonSmooth(dataLeft, dataRight)
            end
          end

          @[AlwaysInline]
          def process(data : Array({{typ.id}})|Slice({{typ.id}})) : Nil
            if @state.smooth?
              processSmooth(data)
            else
              processNonSmooth(data)
            end
          end

          @[AlwaysInline]
          def process(dataLeft : {{typ.id}}, dataRight : {{typ.id}}) : Tuple({{typ.id}}, {{typ.id}})
            left = Slice({{typ.id}}).new(1, dataLeft)
            right = Slice({{typ.id}}).new(1, dataRight)

            if @state.smooth?
              processSmooth(dataLeft, left)
            else
              processNonSmooth(dataLeft, right)
            end

            {left.unsafe_fetch(0), right.unsafe_fetch(0)}
          end
        {% end %}
      {% end %}
    end
  end
end