SIGNAL PROCESSING AND ITS APPLICATIONS

Mariusz Jankowski

University of Southern Maine and Wolfram Research

Abstract

This talk will focus on recent additions and examples of signal processing applications in audio. Recent additions include extended support for Sound and new resampling functionality. Featured applications will include audio filtering, noise removal, voiced/unvoiced detection, frequency shifting, spectral gating, and more.

Application: audio filtering and filter design

Consider the case of a dual-tone multifrequency signal used for telecommunications signaling over analog telephone lines and as used for tone dialing on push-button phones. Each of the keys on the keypad emits a sinusoidal dual-tone from a sum of two sinusoids of two distinct frequencies (in Hz).

DynamicModule[{p,signal},Grid[Join[{{" ",1209,1336,1477}},MapThread[Join,{{{697},{770},{852},{941}},Partition[MapIndexed[Button[First[#2]/.{10->"*",11->0,12->"#"},(EmitSound[Play[Sin[2 Pi #[[1]] t]+Sin[2 Pi #[[2]] t],{t,0,0.2}]]),ImageSize->{40,40}]&,Flatten[Outer[List,{697,770,852,941},{1209,1336,1477}],1]],3]}]]]]

The goal will be to design an analog and digital filter to supress the high frequencies and retain the low.

Application: audio filtering and filter design

Audio visualization - Periodogram

The periodogram is a convenient method of visualizing the frequency content of an audio signal (click on a button to see the spectrum)

DynamicModule[{p,signal},Row[{Grid[Partition[MapIndexed[Button[First[#2]/.{10->"*",11->0,12->"#"},(EmitSound[Play[Sin[2 Pi #[[1]] t]+Sin[2 Pi #[[2]] t]+RandomReal[{-0.2,0.2}],{t,0,0.2}]];
signal=Table[Sin[2 Pi #1[[1]] t]+Sin[2 Pi #1[[2]] t],{t,0,0.2,1/8000.}]+RandomReal[{-0.2,0.2},1601];
p=Periodogram[signal,512,128,BlackmanWindow,SampleRate->8000,PlotRange->All,AspectRatio->2/3,ImageSize->{Automatic,180}]),ImageSize->{40,40}]&,Flatten[Outer[List,{697,770,852,941},{1209,1336,1477}],1]],3]],Spacer[30],Dynamic[p]},Alignment->Top]]

Audio visualization - Spectrogram

The periodogram is a convenient method of visualizing the frequency content of an audio signal (click on a button to see the spectrum)

DynamicModule[{p,signal},Row[{Grid[Partition[MapIndexed[Button[First[#2]/.{10->"*",11->0,12->"#"},(EmitSound[Play[Sin[2 Pi #[[1]] t]+Sin[2 Pi #[[2]] t]+RandomReal[{-0.2,0.2}],{t,0,0.2}]];
signal=Table[Sin[2 Pi #1[[1]] t]+Sin[2 Pi #1[[2]] t],{t,0,0.2,1/8000.}]+RandomReal[{-0.2,0.2},1601];
p=Spectrogram[signal,128,32,BlackmanWindow,SampleRate->8000,PlotRange->All,AspectRatio->2/3,ImageSize->{Automatic,180}]),ImageSize->{40,40}]&,Flatten[Outer[List,{697,770,852,941},{1209,1336,1477}],1]],3]],Spacer[30],Dynamic[p]},Alignment->Top]]

Application: audio filtering and filter design

Here I design two filters that will be able to discriminate between the high and low frequencies of a DTFM signal.

Butterworth filter design

atfB=ButterworthFilterModel[{"Lowpass",{960. 2π,1200. 2π},{1.,30.}},s]

Note that the order of the filter is very high. In practice, this order is too high to be realized as an analog circuit.

SystemsModelOrder[StateSpaceModel[atfB]]

This shows the magnitude spectrum of the filter and the DTFM frequencies.

Plot[atfB[I 2 π f]//Abs,{f,500,2000.},PlotRange->All,Epilog->{Red,Line[{{#,0.},{#,1.}}]&/@{697,770,852,941,1209,1336,1477}},ImageSize->300]

Elliptic filter design

atfE=EllipticFilterModel[{"Lowpass",{960. 2π,1200. 2π},{1.,30.}},s]

Note how the order is significantly lower. This would be easily realizable as an analog circuit.

SystemsModelOrder[StateSpaceModel[atfE]]

This shows a partitioning of the transfer function into first- and secong-order sections (see Wednesday talk by T. Samardzic).

This shows the magnitude spectrum of the filter and the DTFM frequencies.

Plot[atfE[I 2 π f]//Abs,{f,500,2000.},PlotRange->All,Epilog->{Red,Line[{{#,0.},{#,1.}}]&/@{697,770,852,941,1209,1336,1477}},ImageSize->300]

Application: audio filtering and filter design

As we will be filtering a digital version of an example dtfm file we need to convert the analog filter to a digital one. This gives us an infinite impulse response (i.e., a recursive) filter that we will need to apply to the signal using RecurrenceFilter.

Analog-to-digital conversion

Despite the high order of the Butterworth filter, it is preferrable to use it in a digital realization because its phase response is relatively close to linear. This is an important consideration in audio applications.

dtf=ToDiscreteTimeModel[atfB,1/8000.,z]//Chop

This shows the magnitude spectrum on the normalized interval 0≤f≤0.5.

Plot[dtf[(E^(I f 2 π))]//Abs,{f,0.05,0.25},PlotRange->All,Epilog->{Red,Line[{{#/8000.,0},{#/8000.,1}}]&/@{697,770,852,941,1209,1336,1477,1633}},ImageSize->300]

Here is the impulse response of the filter

ListPlot[RecurrenceFilter[dtf,PadRight[{1},100]],PlotRange->All,Filling->0,AspectRatio->1/3,ImageSize->400]

Recursive filtering

This reads an example signal.

snd=;

res=Sound[SampledSoundList[{N@Chop@RecurrenceFilter[dtf,snd[[1,1,1]]]},8000]]

Visualization

Compare the spectrograms of the two signals using a 30 msec window length.

n=Floor[0.03 8000];d=Floor[n/4];
GraphicsRow[{Spectrogram[snd[[1,1,1]],n,d,HammingWindow,SampleRate->8000],Spectrogram[res[[1,1,1]],n,d,HammingWindow,SampleRate->8000]}]

What’s new

ArrayResample - a new function for interpolating sampled data, is one of the more important additions to the signal processing functionality in version 10

This function is critical because of the frequent need to resample audio signals and more generaly, data. There is a multitude of sampling rates used in the audio industry, not including custom rates in use in a variety of databases that store audio data for research purposes.

it is well known that the ideal interpolator is the Sinc function but since it has infinite support it is not used in practice. This illustrates interpolation as a superposition of scaled and shifted sinc functions.

f(x) = (∑_(k=-∞)^∞)f(k) sinc(π (x-k))

GraphicsRow[{Plot[{-0.5 Sinc[π(x-1)],1.5Sinc[π(x-2)],Sinc[π(x-3)]},{x,-2,5},PlotRange->All,Filling->0,Epilog->{Red,{PointSize[Medium],Point[{{1,-0.5},{2,1.5},{3,1.}}]},Line[{{1,0},{1,-0.5}}],Line[{{2,0},{2,1.5}}],Line[{{3,0},{3,1.}}]}],Plot[Total[{-0.5 Sinc[π(x-1)],1.5Sinc[π(x-2)],Sinc[π(x-3)]}],{x,-2,5},PlotRange->All,Filling->0,Epilog->{Red,{PointSize[Medium],Point[{{1,-0.5},{2,1.5},{3,1.}}]},Line[{{1,0},{1,-0.5}}],Line[{{2,0},{2,1.5}}],Line[{{3,0},{3,1.}}]}]},Scaled[0.5],ImageSize->600]

Resampling kernels

In practice resampling functions with finite support are used. Here is a small selection of the resampling kernels available in the ArrayResample function.

GraphicsGrid[Partition[MapThread[Plot[#2[x],{x,-3,3},PlotRange->All,Ticks->{Automatic,False},Exclusions->None,Filling->0,PlotLabel->#1]&,{{"Nearest","Linear","Dodgson","CatmullRom","Schaum","Blackman","Keys","BSplineCubic"},{UnitBox,UnitTriangle,dodgson,catmullrom,schaumquadratic,blackman,keys,schaumcubic}}],4],Spacings->{Automatic,50},Frame->True,ImageSize->900]

ArrayResample

Resampling of sampled data {(Subscript[a, 1]), (Subscript[a, 2]),..., (Subscript[a, n])} is based on the following convolution-like formulation:

f(x) = (∑_(k=1)^n)(Subscript[a, k]) ϕ(x-k)

where function ϕ(x) is a resampling kernel with local support and property ϕ(k)=({{{, ({{1,, k=0}, {0,, k= ±1,±2,±3,...}})}}) .

ArrayResample returns the list {f((Subscript[x, 1])), f((Subscript[x, 2])), ..., f((Subscript[x, m]))} where sample positions (Subscript[x, i]) = (Subscript[x, 0])+((n-1)/m)i  define a uniform sampling of the domain of function f(x) with (Subscript[x, 0]) an arbitrary, user-selected starting position.

data = {0,0,-1,2,3,3,1,4,3,2,1,0,0};
n=13;d=30;
GraphicsGrid[{{Plot[Total@MapIndexed[#1 nearest[x-First@#2]&,data],{x,0,Length[data]+1},PlotRange-> {-1.5,4.5},Filling-> Axis,PlotPoints->100,Axes->False,PlotLabel->"Nearest",
Epilog->{{Red,PointSize[0.02],Point@ MapIndexed[{First@#2,#1}&,data],MapIndexed[Line[{{First@#2,0},{First@#2,#1}}]&,data]},{Blue,PointSize[0.015],Point@Transpose[{Table[0.5(1+n/d)+i n/d,{i,0,d-1}],ArrayResample[N@data,d,Resampling->"Nearest"]}]}}],
Plot[Total@MapIndexed[#1 linear[x-First@#2]&,data],{x,0,Length[data]+1},PlotRange-> {-1.5,4.5},Filling-> Axis,PlotPoints->100,Axes->False,PlotLabel->"Linear",
Epilog->{{Red,PointSize[0.02],Point@ MapIndexed[{First@#2,#1}&,data],MapIndexed[Line[{{First@#2,0},{First@#2,#1}}]&,data]},{Blue,PointSize[0.015],Point@Transpose[{Table[0.5(1+n/d)+i n/d,{i,0,d-1}],ArrayResample[N@data,d,Resampling->"Linear"]}]}}],Plot[Total@MapIndexed[#1 dodgson[x-First@#2]&,data],{x,0,Length[data]+1},PlotRange-> {-1.5,4.5},Filling-> Axis,PlotPoints->100,Axes->False,PlotLabel->"Dodgson",
Epilog->{{Red,PointSize[0.02],Point@ MapIndexed[{First@#2,#1}&,data],MapIndexed[Line[{{First@#2,0},{First@#2,#1}}]&,data]},{Blue,PointSize[0.015],Point@Transpose[{Table[0.5(1+n/d)+i n/d,{i,0,d-1}],ArrayResample[N@data,d,Resampling->"Dodgson"]}]}}],
Plot[Total@MapIndexed[#1 catmullrom[x-First@#2]&,data],{x,0,Length[data]+1},PlotRange-> {-1.5,4.5},Filling-> Axis,PlotPoints->100,Axes->False,PlotLabel->"CatmullRom",
Epilog->{{Red,PointSize[0.02],Point@ MapIndexed[{First@#2,#1}&,data],MapIndexed[Line[{{First@#2,0},{First@#2,#1}}]&,data]},{Blue,PointSize[0.015],Point@Transpose[{Table[0.5(1+n/d)+i n/d,{i,0,d-1}],ArrayResample[N@data,d,Resampling->"CatmullRom"]}]}}]},
{Plot[Total@MapIndexed[#1 schaumcubic[x-First@#2]&,data],{x,0,Length[data]+1},PlotRange-> {-1.5,4.5},Filling-> Axis,PlotPoints->100,Axes->False,PlotLabel->"Schaum",
Epilog->{{Red,PointSize[0.02],Point@ MapIndexed[{First@#2,#1}&,data],MapIndexed[Line[{{First@#2,0},{First@#2,#1}}]&,data]},{Blue,PointSize[0.015],Point@Transpose[{Table[0.5(1+n/d)+i n/d,{i,0,d-1}],ArrayResample[N@data,d,Resampling->"SchaumCubic"]}]}}],Plot[Total@MapIndexed[#1 blackman[x-First@#2]&,data],{x,0,Length[data]+1},PlotRange-> {-1.5,4.5},Filling-> Axis,PlotPoints->100,Axes->False,PlotLabel->"Blackman",
Epilog->{{Red,PointSize[0.02],Point@ MapIndexed[{First@#2,#1}&,data],MapIndexed[Line[{{First@#2,0},{First@#2,#1}}]&,data]},{Blue,PointSize[0.015],Point@Transpose[{Table[0.5(1+n/d)+i n/d,{i,0,d-1}],ArrayResample[N@data,d,Resampling->"Blackman"]}]}}],Plot[Total@MapIndexed[#1 keys[x-First@#2]&,data],{x,0,Length[data]+1},PlotRange-> {-1.5,4.5},Filling-> Axis,PlotPoints->100,Axes->False,PlotLabel->"Keys",
Epilog->{{Red,PointSize[0.02],Point@ MapIndexed[{First@#2,#1}&,data],MapIndexed[Line[{{First@#2,0},{First@#2,#1}}]&,data]},{Blue,PointSize[0.015],Point@Transpose[{Table[0.5(1+n/d)+i n/d,{i,0,d-1}],ArrayResample[N@data,d,Resampling->"Keys"]}]}}],Plot[Total@MapIndexed[#1 bsplinecubic[x-First@#2]&,filteredDataCubic],{x,0,Length[data]+1},PlotRange-> {-1.5,4.5},Filling-> Axis,PlotPoints->100,Axes->False,PlotLabel->"BsplineCubic",
Epilog->{{Red,PointSize[0.02],Point@ MapIndexed[{First@#2,#1}&,data],MapIndexed[Line[{{First@#2,0},{First@#2,#1}}]&,data]},{Blue,PointSize[0.015],Point@Transpose[{Table[0.5(1+n/d)+i n/d,{i,0,d-1}],ArrayResample[N@data,d,Resampling->"BSplineCubic"]}]}}]}},Frame->True,ImageSize->900]

ArrayResample

Interpolation with cubic Bspline kernel.

data = {0,0,-1,2,3,3,1,4,3,2,1,0,0};
n=13;d=19;
Plot[Total@MapIndexed[#1 bsplinecubic[x-First@#2]&,filteredDataCubic],{x,0,Length[data]+1},PlotRange-> {-1.5,4.5},Filling-> Axis,PlotPoints->100,
Epilog->{{Red,PointSize[0.02],Point@ MapIndexed[{First@#2,#1}&,data],MapIndexed[Line[{{First@#2,0},{First@#2,#1}}]&,data]},{Blue,PointSize[0.015],Point@Transpose[{Table[0.5(1+n/d)+i n/d,{i,0,d-1}],ArrayResample[N@data,19,Resampling->"BSplineCubic"]}]}}]

Vary the number of sample points and starting positions.

DynamicModule[
{data,n,x,filteredData},
data={0,0,-1,2,3,3,1,4,3,2,1,0,0};
n=Length[data];
Manipulate[
Plot[Total@MapIndexed[#1 bsplinecubic[(x-First@#2)]&,filteredData],{x,If[pos<0,pos,0],Length[data]+1},
PlotRange-> {-1.5,4.5},Filling->0,PlotPoints->100,AxesOrigin->{0,0},
Epilog->{{Red,PointSize[0.02],Point@ MapIndexed[{First@#2,#1}&,data]},{Blue,PointSize[0.015],Point[Thread[{Table[pos+i n/d,{i,0,d-1}],ArrayResample[data,d,pos,Resampling->"BSplineCubic"]}]]}}],
{{d,13,"number of samples"},13,36,1},{{pos,0.5(1+n/d),"starting position"},-2,4},
Initialization:>(filteredData={-0.17131897704272686`,0.6393711266241827`,-2.3861655294540043`,2.9052909911918356`,2.7650015646866626`,4.034702750061514`,-0.90381256493272`,5.580547509669365`,2.5816225262552592`,2.0929623853095953`,1.0465279325063581`,-0.2790741153350289`,0.06976852883375725`};)]
]

Application: audio resampling

Upsample audio signal using different interpolating kernels.

snd=Import["C:\\Users\\Mariusz\\Documents\\Data\\audio\\rampal[1].wav"];
data=snd[[1,1,1]][[Floor[2.5 srate];;Floor[2.5 srate]+30799]];srate=snd[[1,2]];
{ListPlay[data,PlayRange->{-1,1},SampleRate->srate],ListPlay[ArrayResample[data,4Length[data],1,Resampling->"Linear"],PlayRange->{-1,1},SampleRate->4 srate],ListPlay[ArrayResample[data,4Length[data],1,Resampling->"Keys"],PlayRange->{-1,1},SampleRate->4 srate]}

Downsample audio signal using different interpolating kernels.

snd=Import["C:\\Users\\Mariusz\\Documents\\Data\\audio\\cwarb[1].wav"];
data=snd[[1,1,1]][[Floor[0.5 srate];;Floor[0.5 srate]+30799]];srate=snd[[1,2]];
{ListPlay[data,SampleRate->srate],ListPlay[ArrayResample[data,Length[data]/4,1,Resampling->"Linear"],PlayRange->{-1,1},SampleRate->srate/4],ListPlay[ArrayResample[data,Length[data]/4,1,Resampling->"Keys"],PlayRange->{-1,1},SampleRate->srate/4]}

ArrayResample and upsampling/convolution

An alternative interpretation of interpolation as a two step operation of upsampling and convolution.

data={0,0,-1,2,3,3,1,4,3,2,1,0,0};

Upsample[{0,0,-1,2,3,3,1,4,3,2,1,0,0},2]

ListConvolve[{(1/2),1,(1/2)},Upsample[{0,0,-1,2,3,3,1,4,3,2,1,0,0},2],2,0]

Compare with ArrayResample

ArrayResample[{0,0,-1,2,3,3,1,4,3,2,1,0,0},26,1,Resampling->"Linear"]

Nyquist filters

While the multitude of built-in resampling kernels provide adequate flexibility for most resampling tasks, on occasion it may be necessary to have greater control over the frequency response of an interpolating filter. Here I show how to design Nyquist filters. A Nyquist filter of order k has a cutoff frequency at (π/k).

We can estimate the needed filter length using Kaiser’s formula. Assuming we want stopband attenuation of δ=0.01 (or -40 dB) and a transition band width of Δω=π/8 rad/sec we get

((20 Log10[100]-13)/(2.324 (π/8)))

This gives the desired filter

fir=N@Chop@LeastSquaresFilterKernel[{"Lowpass",π/4},31] Array[HammingWindow,31,{-0.5,0.5}]

This resamples the audio using the designed filter.

{ListPlay[data,SampleRate->srate],ListPlay[ListConvolve[fir,Upsample[data,4],16,0.],SampleRate->4 srate]}

Compare frequency response of three filters/kernels: linear, Lanczos and the Nyquist filter fir.

linear=((#)/(Total[#]))&[Table[UnitTriangle[x],{x,-1,1,0.25}]];
keys=((#)/(Total[#]))&[Table[keys[x],{x,-2,2,0.25}]];
Plot[Evaluate[20Log10[Abs[ListFourierSequenceTransform[#,ω]&/@{linear,keys,fir}]]],{ω,0,π},GridLines->Automatic,PlotLegends->{"linear","Lanczos","Nyquist"}]

Application: frequency shifting

Frequency shift theorem.

Given a periodic sequence x[n] (with period N) and its discrete-time Fourier series X[k]  we have

x[n]  ·(e^(j ((2 π)/N)n (Subscript[k, 0])))  (Overscript[⟷, DFT])  X[k-(Subscript[k, 0])]

Here is an example for (Subscript[k, 0])=2

Fourier[Table[Cos[((2π)/9)n],{n,0,8}],FourierParameters->{-1,-1}]//Chop

Fourier[Table[Cos[((2π)/9)n] Exp[I ((2π)/9)n] ,{n,0,8}],FourierParameters->{-1,-1}]//Chop

Observe the Fourier spectrum as the shift value (Subscript[k, 0]) is varied. Clearly, the “shift” is in effect a rotation within the list of Fourier coefficients returned by function Fourier.

Manipulate[ListPlot[Abs[Fourier[Table[(Cos[((2π)/63)7n]+(1/2)Sin[((2π)/64)13n]) Exp[I ((2π)/64)(Subscript[k, 0]) n] ,{n,0,63}],FourierParameters->{-1,-1}]],Filling->0,AxesOrigin->{0,0},PlotRange->All,PlotStyle->PointSize[Medium],AspectRatio->1/2],{(Subscript[k, 0]),0,19,1}]

Application: frequency shifting

This defines an example audion signal.

digits=;

Typical frame duration for speech processing is10-40 msec. Here I use 20 msec to calculate the spectorgram.

data=digits[[1,1,1]];srate=digits[[1,2]];
n=If[OddQ[#],#,#+1]&[Floor[srate 0.02]];d=Floor[n/3];
sa=SpectrogramArray[data,n,d];frames=First@Dimensions[sa];
(*Print[{n,d,frames,Round[srate/n]}];*)
Spectrogram[data,n,d,SampleRate->srate]

Up-shift frequencies by 550 Hz.

result=Flatten[Take[InverseFourier[fshift[#,Floor[550 n/srate]],FourierParameters->{1,-1}],d]&/@sa];
ListPlay[result,PlayRange->{-1,1},SampleRate->srate]

Down-shift frequencies by 250 Hz.

result=Flatten[Take[InverseFourier[fshift[#,-Floor[250 n/srate]],FourierParameters->{1,-1}],d]&/@sa];
ListPlay[result,PlayRange->{-1,1},SampleRate->srate]

Application: voiced/unvoiced detection

The problem of locating the beginning and end of a speech utterance in a background of noise (or other acoustic signals) is important in many areas of speech processing. In particular, in automatic speech recognition, it is essential to locate the time intervals when speech is present so as to avoid trying to classify background sounds as speech inputs. A scheme for locating the beginning and end of a speech signal can also be used to eliminate signiﬁcant computation by making it possible to process only the parts of an input signal that correspond to speech.

digits=;

A typical frame duration used in speech processing is 10-40 msec. Here I use 40 msec.

{n=Floor[0.04 44100],d=Floor[n/4]}

Partition the speech signal into overlapping frames.

data=digits[[1,1,1]];
win=Array[HammingWindow,n,{-0.5,0.5}];
arrwin=Map[win #&,Partition[data,n,d]];

Short-time energy

The total energy of voiced speech is considerably larger than that of unvoiced speech. Therefore, voiced speech should be characterized by relatively high energy while unvoiced speech will have relatively low energy. This shows the short-time energy of the example audio signal.

ListLinePlot[energy=Table[Total@((arrwin[[r]])^2),{r,1,Length[arrwin]}],PlotRange->All,ImageSize->Small]

This applies nonlinear smoothing to the energy signal and calculates a threshold based on speech processing heuristics to detect voiced/unvoiced segments.

er=10Log10[#]-10Log10[Max[#]]&[energy];
m=Mean[er[[1;;20]]]; s=StandardDeviation[er[[1;;20]]];te=Min[-30,m+3 s];
ListLinePlot[MedianFilter[er,3],PlotRange->All,Epilog->{Red,Line[{{1,te},{Length[arrwin],te}}]},ImageSize->300]

Short-time zero-crossings

The model for speech production suggests that for unvoiced speech, most of the energy is found at higher frequencies. Since high frequency spectral content implies high zero-crossing rates, and low frequency spectral content implies low zero-crossing rates, there is a strong correlation between zero-crossing rate and energy distribution with frequency. Therefore, voiced speech should be characterized by relatively low zero-crossing rate, while unvoiced speech will have a relatively high zero-crossing rate.

zero=N[(d/(2n))Total[CrossingDetect[#,CornerNeighbors->None]]&/@arrwin];
m=Mean[zero[[1;;20]]];s=StandardDeviation[zero[[1;;20]]]; tz=Max[35,m+3 s];
ListLinePlot[MedianFilter[zero,3],PlotRange->All,Epilog->{Red,Line[{{1,tz},{Length[arrwin],tz}}]},ImageSize->300]

Combine the two measures

emask=If[#<te,0,1]&/@MedianFilter[er,3];
zmask=If[#>tz,1,0]&/@MedianFilter[zero,3];
ListLinePlot[{zmask,emask},ImageSize->300]

Application: spectrogram modulation

In this example I modify an audio signal by modulating it in the frequency domain with selected images. This shows a sample of audio and its spectrogram.

snd=;
data=snd[[1,1,1]];srate=snd[[1,2]];
n=If[OddQ[#],#,#+1]&[Floor[0.04 srate]];d=Floor[n/4];
sa=SpectrogramArray[data,n,d];{h,w}=Dimensions[sa];
Row[{snd,Spectrogram[data,n,d,HannWindow,SampleRate->srate,ImageSize->400]},Spacer[50]]

The basic idea is to multiply the magnitude of the short-time Fourier transform (i.e., the spectrogram array) by a grayscale image. We need to impose even symmetry on the modulating image to ensure that the returned signal is real.

“cloudy” Rampal

This creates an image with even symmetry.

i=ImageAdjust@ImageAssemble[{{#,ImageReflect[ImageTake[#,All,{2,-1}],Left]}}]&[GradientFilter[ColorConvert[ImageTake[,{1,h},{1,((w-1)/2)+1}],"GrayLevel"],3]]

Perform the modulation and invert the spectrogram to obtain the resulting sound.

modulatedSound=Re@Flatten[Take[InverseFourier[#,FourierParameters->{1,-1}],d]&/@(Abs[sa] ImageData[i] Exp[I Arg[sa]])];
Row[{ListPlay[modulatedSound,SampleRate->srate],Spectrogram[modulatedSound,n,d,SampleRate->srate,ImageSize->500]},Spacer[15]]

“Athenian” Rampal

i=ImageAdjust@ImageAssemble[{{#,ImageReflect[ImageTake[#,All,{2,-1}],Left]}}]&[GradientFilter[ImageResize[,{((w-1)/2)+1,h}],3]]

modulatedSound=Re@Flatten[Take[InverseFourier[#,FourierParameters->{1,-1}],d]&/@(Abs[sa] ImageData[i,Real] Exp[I Arg[sa]])];
Row[{ListPlay[modulatedSound,SampleRate->srate],Spectrogram[modulatedSound,n,d,SampleRate->srate,ImageSize->500]},Spacer[15]]

Application: the “sound” of graphics and images

octopus

oct=ImageAssemble[{#,ImageReflect[ImageTake[#,All,{2,-1}],Left]}]&[ImageAdjust[ImageReflect[ImageRotate[GaussianFilter[ColorConvert[ImageResize[,Scaled[0.25]],GrayLevel],11],Right],Left]]]

{w,h}=ImageDimensions[oct];
sound=Re@Flatten[Take[InverseFourier[#,FourierParameters->{1,-1}],Floor[w/4]]&/@ImageData[oct]];
Row[{ListPlay[sound,SampleRate->22050],Spectrogram[sound,w,Floor[w/4],SampleRate->22050,ImageSize->400]},Spacer[50]]

Thank you

Questions?

Initialize

Define resampling kernels

ideal[x_]:=Sinc[π x];
nearest[x_]:=UnitBox[x];
linear[x_]:=UnitTriangle[x];
dodgson[x_]:=Piecewise[{{3/2+(5*x)/2+x^2,-3/2<x<-1/2},{1/2,x==-1/2},{1-2*x^2,-1/2<x<1/2},{1/2,x==1/2},{3/2-(5*x)/2+x^2,1/2<x<3/2}},0];
catmullrom[x_]:=Piecewise[{{2+4*x+(5*x^2)/2+x^3/2,-2<x<-1},{1-(5*x^2)/2-(3*x^3)/2,-1<x<0},{1,x==0},{1-(5*x^2)/2+(3*x^3)/2,0<x<1},{2-4*x+(5*x^2)/2-x^3/2,1<x<2}},0];
german[x_]:=(1/144)*Piecewise[{{(-(3+x))*(4+x)^2,-4<x≤-3},{-48-4*x+40*x^2+21*x^3+3*x^4,-3<x≤-2},{312+580*x+306*x^2+29*x^3-9*x^4,-2<x≤-1},{144-335*x^2-185*x^3+6*x^4,-1<x≤0},{144-335*x^2+185*x^3+6*x^4,0<x≤1},{312-580*x+306*x^2-29*x^3-9*x^4,1<x≤2},{-48+4*x+40*x^2-21*x^3+3*x^4,2<x≤3},{(-4+x)^2*(-3+x),3<x<4}},0];
schaumquadratic[x_]:=Piecewise[{{-1/16,x==-3/2},{1+(3*x)/2+x^2/2,-3/2<x<-1/2},{9/16,x==-1/2},{1-x^2,-1/2<x<1/2},{9/16,x==1/2},{1-(3*x)/2+x^2/2,1/2<x<3/2},{-1/16,x==3/2}},0];
schaumcubic[x_]:=Piecewise[{{1+(11*x)/6+x^2+x^3/6,Inequality[-2,Less,x,LessEqual,-1]},{1+x/2-x^2-x^3/2,Inequality[-1,Less,x,LessEqual,0]},{1-x/2-x^2+x^3/2,Inequality[0,Less,x,LessEqual,1]},{1-(11*x)/6+x^2-x^3/6,1<x<2}},0];
keys[x_]:=Piecewise[{{0,x==-2},{2+4*x+(5*x^2)/2+x^3/2,-2<x<-1},{0,x==-1},{1-(5*x^2)/2-(3*x^3)/2,-1<x<0},{1,x==0},{1-(5*x^2)/2+(3*x^3)/2,0<x<1},{0,x==1},{2-4*x+(5*x^2)/2-x^3/2,1<x<2}},0];
blackman[x_]:=BlackmanWindow[x/8]*Sinc[Pi*x];
hamming[x_]:=HammingWindow[x/8]*Sinc[Pi*x];
lanczos[x_] := LanczosWindow[x/6] Sinc[π x];
bsplinequadratic[x_] := Piecewise[{{(3+2*x)^2,-3/2<x≤-1/2},{6-8*x^2,-1/2<x≤1/2},{(3-2*x)^2,1/2<x<3/2}},0]/8;
bsplinecubic[x_] := Piecewise[{{(2+x)^3/6,-2<x≤-1},{(4-6*x^2-3*x^3)/6,-1<x≤0},{(4-6*x^2+3*x^3)/6,0<x≤1},{-(-2+x)^3/6,1<x<2}},0];

inverseQuadraticBsplineFilter[data_]:= Reverse@RecurrenceFilter[(TransferFunctionModelFull[(8/(1.0+(3+2 (Sqrt[2])) z)), 1, , ]),Reverse@RecurrenceFilter[(TransferFunctionModel[(1/(3.0-2(Sqrt[2])+z)), 1]),data]];

filteredDataQuadratic=inverseQuadraticBsplineFilter[{0,0,-1,2,3,3,1,4,3,2,1,0,0}]

inverseCubicBsplineFilter[data_]:= Reverse@RecurrenceFilter[(TransferFunctionModel[(6/(1.0+(2+(Sqrt[3])) z)), 1]),Reverse@RecurrenceFilter[(TransferFunctionModel[(1/(2.0-(Sqrt[3])+z)), 1]),data]];

filteredDataCubic=inverseCubicBsplineFilter[{0,0,-1,2,3,3,1,4,3,2,1,0,0}]

fshift[x_,k_]:=Module[{n,y,z,dc},
n=Length[x];
If[EvenQ[n],
Return[$Failed],
y=x;dc=y[[1]];
If[k>0,
(y[[-k;;-1]]=0.;
y=Fourier[InverseFourier[y,FourierParameters->{1,-1}]Exp[I 2 π/n k Range[0.,n-1.]],FourierParameters->{1,-1}];
y[[1;;k+1]]=0.;
y[[1]]=dc;),
(y=Fourier[InverseFourier[y,FourierParameters->{1,-1}]Exp[I 2 π/n k Range[0.,n-1.]],FourierParameters->{1,-1}];
y[[1]]=dc;)
];
Join[#,Conjugate[Reverse[Rest[#]]]]&[y[[1;;(n+1)/2]]]
]
]

do-loop code for frequency shift

(*
shift=Floor[550 n/srate];
n2=(n-1)/2;new=sa;
Do[tmp=new[[j,1;;n2]];
Do[tmp[[n2-i+1]]=tmp[[n2-i+1-shift]],{i,1,n2-shift-1}];
tmp[[2;;shift]]=0.;
new[[j]]=Join[tmp,{sa[[j,n2]]},Conjugate[Reverse[Rest[tmp]]]],{j,1,frames}];
result=Re@Flatten[Take[#,d]&/@(InverseFourier[#,FourierParameters->{1,-1}]&/@new)];
ListPlay[result,PlayRange->{-1,1},SampleRate->srate]
*)