• Welcome to ASR. There are many reviews of audio hardware and expert members to help answer your questions. Click here to have your audio equipment measured for free!

Some elementary Matlab programs to test your system.

Thanks for the detailed reply. I've not looked much into the theory behind Heyser's time delay spectrometry and related techniques. Looks like I have some reading to do :).

What was the intended subject?
I believe the poster was pointing out that REW is capable enough to replace a number of special-purpose commercial measurement/analysis programs. Arguably not directly relevant to this thread, though it did serve as a prompt to explaining why someone might want to write their own tooling (which was unclear to some readers). I have no idea where the "marketing" allegations came from...
 
  • Like
Reactions: j_j
Test your MATLAB system with basic programs like displaying "Hello, World", performing simple arithmetic operations, and plotting a sine wave. These check core functions: output, computation, and graphics rendering.
 
Test your MATLAB system with basic programs like displaying "Hello, World", performing simple arithmetic operations, and plotting a sine wave. These check core functions: output, computation, and graphics rendering.

"Hello World"?

sprintf('Hello World \n')
2+2
plot(sin(2*pi*(0:512)/128));

Enjoy. Note no ';' in the first two lines.

Now, MATLAB isn't "mine" so I don't understand your use of "your", but do you actually know what MATLAB is?
 
  • Like
Reactions: NTK
I finally write a quick script to analyze the horrible buzz tone.

first, 3 jpg's
nodis.jpg


interch.jpg
dist2nd.jpg

And the script. Note, file names may mislead, except that the very first pair are undistorted l and r.

clc
clear
close all

len=32768;

%run the nasty buzztone through your system
%loop it many times, without windowing,
%and it will be continuous.
%after everything is in steady state, grab exactly 32768 samples.
%or grab more and then truncate to exactly 32768 contiguous samples.
%put that into ergh.wav (or set your own filename, whatever)

[x, fs]=audioread('ergh.wav');

xt=fft(x,len);

xta=abs(xt(1:(len/2+1),:));

stmax=max(max(xta)); % assuming two or more channels
xta=xta/stmax; %normalize
xta=max(xta,.000001); % no log(zero) please

subplot(2,1,1)
semilogy(xta(:,1))
subplot(2,1,2);
semilogy(xta(:,2))

Anyone want to guess what the second and third plots are? Yes, they both result from 2nd order nonlinearities. The reason for this buzz tone is that it creates those lines, inexorably, between the tones in the original.
 
Trying to follow. Have no idea what kinds of nonlinearities the second and third plots are showing :facepalm:

Here is my attempt to port J J's code to Julia, since I am not much of a MATLAB user and I don't have the DSP/audio toolkits. The code to generate and plot the buzz tone.
Code:
using Plots;
using Primes;
using WAV;
using FFTW;
using Random;
pythonplot(size=(1000, 500));    # Use the pythonplot backend for plotting

lfreq = prevprime(250);    # Want the two prime numbers closest to 250, one above and one below
rfreq = nextprime(250);
fs = 48000;

len = 2^15;
rng = Xoshiro(1234);    # Use the Julia default Xoshiro256++ random number generator, with given seed

ilf = round(Int32, len/fs*lfreq - 1);
irf = round(Int32, len/fs*rfreq - 1);
istride = round(Int32, len/fs*500);    # Yes, this partially obviates the use of primes above, but useful as we shall see

xt = zeros(ComplexF64, (len, 2));
xt[ilf:istride:(len÷2), 1] .= 1;
xt[irf:istride:(len÷2), 2] .= 1;
for ii in 2:len÷2
    global t = rand(rng)*2π;
    global xt[ii, 1:2] *= cos(t) + 1im*sin(t);   # Reduce peak to rms ratio more or less
end;

xt[len:-1:(len÷2 + 2), 1:2] = conj.(xt[2:(len÷2), 1:2]);
x = real(ifft(xt));

t = maximum(maximum(abs.(x)));
x *= 0.9/t;

plt = plot(x);

for ii in 1:6    # Increase length by catenation (hence the FFT use for glitchfree outcome)
    global x = [x; x];
end;

filename = joinpath(homedir(), "JJ_buzz_tone.wav");
wavwrite(x, filename, Fs=fs);

savefig(plt, joinpath(homedir(), "Pictures", "buzz_tone_waveforms.png"));
plot(plt)
Here is the plot of the buzz tone waveform.
buzz_tone_waveforms.png


Code to read back the WAV file and plot the spectra.
Code:
using Plots;
using WAV;
using FFTW;
using Printf;
pythonplot(size=(1000, 500));    # Use the pythonplot backend for plotting

abs_clipped(x) = min( max(abs(x), 1e-12), 1e12 );    # Return the magitude of a complex number, clip to 1e-12 ≤ |x| ≤ 1e12

filename = joinpath(homedir(), "JJ_buzz_tone.wav");
x, fs = wavread(filename);
nsamples = size(x)[1];

nfftpts = 2^15;
hx = rfft(x[1:nfftpts, :])/(2π*sqrt(nfftpts))

xticks = 0:2000:18000;

p1 = plot(abs_clipped.(hx[:, 1]), yscale=:log10, xlims=(0, xticks[end]), ylims=(1e-6, 2), label="Ch A");
xticks!( p1, ( xticks, [ @sprintf("%.0f", xt) for xt in xticks ] ) );

p2 = plot(abs_clipped.(hx[:, 2]), yscale=:log10, xlims=(0, xticks[end]), ylims=(1e-6, 2), label="Ch B");
xticks!( p2, ( xticks, [ @sprintf("%.0f", xt) for xt in xticks ] ) );

plt = plot(p1, p2, layout=(2, 1));
savefig(plt, joinpath(homedir(), "Pictures", "buzz_tone_spectra.png"));
plot(plt)
buzz_tone_spectra.png

Eagerly awaiting J J to reveal the secrets, or provide some hints :)
 
  • Like
Reactions: j_j
The second of the three plots was a cross-channel multiply with a gain of .1 (i.e. L = L + .1 L*R) The third was simple 2nd order (L + .1*L*L). both per-sample operations.
 
  • Like
Reactions: NTK
The pgen.m script you've provided is designed to generate a broadband noise signal with a uniform spectral distribution, resulting in a 10-second audio file named whooosh.wav. This type of signal is particularly useful for audio system testing, psychoacoustic experiments, or as a creative sound design element.


Key Features of the Script:


  • Impulse Initialization: The script starts by creating an impulse at the center of a zero-initialized array, ensuring a wide frequency content upon transformation.
  • Phase Randomization: A random phase spectrum is generated and smoothed using a low-pass Butterworth filter. This process ensures a natural-sounding variation in the noise while maintaining a flat frequency response.
  • Spectrum Symmetry: To guarantee a real-valued time-domain signal after the inverse FFT, the script constructs a conjugate-symmetric frequency spectrum.
  • Normalization: The resulting signal is normalized to prevent clipping and to utilize the full dynamic range of the audio format.
  • Output: The final audio is saved as a 24-bit WAV file, providing high-resolution audio suitable for professional applications.

Almost reads like AI but it's right, so it's not AI. You want a job? :)


You might mention as well it makes deconvolution via FFT extremely easy.
 
This creates a bit longer than 10 second clip of weird broadband noise. It's value is that it is energetic, and has no zeros or near zeros in its spectrum.
It creates a file called "whooosh.wave".

You may have to swap to wavwrite on octave. The rest should "just work" in octave, although I lack a modern octave any longer, having actual matlab.

---file pgen.m---
clc
clear all
close all

fs=48000;

blen=2^ceil(log2(10*fs)); %next step up from 10 seconds.

x(1:blen,1)=0;
x(round(blen/2),1)=1;

xt=fft(x);

%phi only needs to be 2:blen/2 memory is cheap
phi=rand(blen,1)*2*pi;% many will be unused, so?
phi=phi-sum(phi)/blen; % make DC value zero
phi=phi *5.6 ; %just fiddle for the right spread


[bb,aa]=butter(4, .1); %the fractional lp frequency affects spread too

phi=filter(bb,aa,phi); %limit speed of phase change.

phi(1)=0; % to prevent complex DC part

for ii=2:blen/2
xt(ii)=xt(ii)*( cos(phi(ii)) + i*sin(phi(ii)));
end
xt(blen:-1: (blen/2+2))=conj(xt(2: (blen/2))); % conjugate spectrum for real only output.
%notice 1 and blen/2+1 are unmodified and purely real, not complex

xx=ifft(xt);
t=max(abs(xx));
xx=xx*.9/t;

audiowrite('whooosh.wav',xx,fs,'BitsPerSample', 24);

plot(xx)
figure
plot(abs(fft(xx,2*blen)))
---end pgen.m---

Free for non-commercial use.
Nice try FBI.
 
Maybe 5 years ago, I created a Matlab script to reproduce my personal tinnitus experience. Tuned it to blend right in with the tones, mostly 8-16 kHz with adjacent tones with annoying intermodulation and beats. Tried to play it last year and didn't hear anything. While tinkering with my laptop's audio output settings, my wife and kids started furiously texting me to stop it with the high pitched noises. For me, all those sounds... basically gone, or at least indistinguishable from my ears' noise floor.
 
Well tinnitus usually precedes actual loss. I understand your frustration, I'm also old.
 
165 USD / year ? Kinda expensive, also given the product is not commercial friendly, nor
 
Fun little program to generate some interesting bass.

clc
clear all
close all

fs=48000;

tfreq=50;
len=2^(ceil(log2(4*fs)));

ifreq=round(tfreq*len/fs);



xt(1:len,1:2)=0;

xt(ifreq,1)=1;
xt(ifreq+1,2)=1;

xt(len:-1: (len/2+2),1:2)=conj(xt(2: (len/2),1:2));

x=ifft(xt);


t=max(max(abs(x)))
x=x/t*.9;
plot(x)

audiowrite('weirdbass.wav',x,fs, 'BitsPerSample',16)
 
Fun little program to generate some interesting bass.

clc
clear all
close all

fs=48000;

tfreq=50;
len=2^(ceil(log2(4*fs)));

ifreq=round(tfreq*len/fs);



xt(1:len,1:2)=0;

xt(ifreq,1)=1;
xt(ifreq+1,2)=1;

xt(len:-1: (len/2+2),1:2)=conj(xt(2: (len/2),1:2));

x=ifft(xt);


t=max(max(abs(x)))
x=x/t*.9;
plot(x)

audiowrite('weirdbass.wav',x,fs, 'BitsPerSample',16)

Could someone kindly run this and share the weirdbass.wav file?
 
Back
Top Bottom