BUG: freqz does not work on multi-dimensional array, despite what documentation claims

Describe your issue.

freqz documentation seems to imply that it supports multidimensional inputs as long as the first dimension is the filter coeffients:

b: array_like

Numerator of a linear filter. If b has dimension greater than 1, it is assumed that the coefficients are stored in the first dimension, and b.shape[1:], a.shape[1:], and the shape of the frequencies array must be compatible for broadcasting.

However, freqz fails when given a 2-dimensional array.

Reproducing Code Example

import scipy
import numpy as np

scipy.signal.freqz(np.zeros((128, 10)))

Error message

Traceback (most recent call last):
  File "<string>", line 1, in <module>
  File "/home/sami/.virtualenvs/tmp/lib/python3.10/site-packages/scipy/signal/_filter_design.py", line 474, in freqz
    h = (npp_polyval(zm1, b, tensor=False) /
  File "/home/sami/.virtualenvs/tmp/lib/python3.10/site-packages/numpy/polynomial/polynomial.py", line 754, in polyval
    c0 = c[-1] + x*0
ValueError: operands could not be broadcast together with shapes (10,) (512,)

SciPy/NumPy/Python version information

1.9.3 1.23.4 sys.version_info(major=3, minor=10, micro=6, releaselevel='final', serial=0)