What does Smoothness Prior mean in the context of image processing? I can quote from an article:

  "The first category of methods, known as diffusion-based inpainting,
  introduces smoothness priors via parametric models or partial
  differential equations (PDE) to propagate or diffuse local structures
  from the exterior to the interior of the hole."

What is Smoothness Prior?

What does this definition mean?

A prior is something that you know a priori (beforehand) about how the system you're looking at works.

A smoothness prior means that you know that the image you're looking at will be smooth in a particular way.

That paper goes on to say:

  The term diffusion comes from the idea of propagating local information with smoothness constraints, by analogy with physical phenomena like heat propagation in physical structures. These phenomena can be formalized with PDEs, and diffusion is therefore performed using PDE-based regularization.

So in this case the smoothness that you know about beforehand is enforced using a partial differential equation.


How to Reverse Color Quantization? Are there any algorithms that attempt to perform the inverse of color quantization on an image? In other words, is there a smart way to increase the bit depth of image?

The obvious answer would be to simply stretch the values of all pixels to fit the new range; however, that would leave gaps in the histogram, and the resulting image would still have the same number of different values. Are there any algorithms that, for example, take neighboring pixels into consideration in order to determine a better estimation for each pixel?

Are there any algorithms that, for example, take neighboring pixels into consideration in order to determine a better estimation for each pixel?

That would essentially be a low-pass filter.
So, yes, that exists, and is commonly used.
You can see that very nicely if you take e.g. an old computer graphics sprite, and scale it as a high-color-resolution image:
original (32 x 32 px)

scaled without interpolation filter (512 px)

scaled by a simple linear interpolation filter (512 px)

While the simply scaled image has as much discrete colors as the original, the linearly interpolated one has a lot more colors, and the histogram looks much more "continuous".
But:
There's no free lunch, usually. You can't just add color information back! (not being able to add back lost info is one of the fundamental truths of Information Theory, by the way). You have to sacrifice spatial resolution. Essentially, it's the same math that underlies Heisenberg's Uncertainty Principle. You can either have perfect info state of something (the color of a pixel, the impulse of a electron), or its position (sharp edges in an image, knowledge of the position of an electron).  See this example:
original photo

original (512 x 512px) reduced to 16 colors

filtered to achieve a much more continuous color histogram, with a radius of 10px taking into consideration


Detect non rectangular bounding box for objects I have a short video where a TV is shown. I need to detect the TV there (let's assume it's all pure black or it's covered with a green screen), and overlay a video on top of it on the same area. But the problem is when the angle is not straight, it is oriented and the bounding boz is not a rectangle. So I guess I need to detect the angle and its scale to overlay the content correctly and nicely.
Is there any ways to do it with AI or computer vision? We can detect the TV with object detection, but I don't know how to compute the tilt and angle for best fit when it's not a retangle.

Let's label your corners like this:

1-------2
|       |
3-------4

For your destination.  Let, $P_1, P_2, P_3, \text{ and } P_4 $ be vectors representing the corners in pixel coordinates.
For your source, Let $H$ and $W$ represent the height and width, and $(x,y)$ be a point in that rectangle you want to project.
You can calculate your coefficients like this:
$$
\begin{aligned}
C_1 &= \left(1 - \frac{x}{W} \right)\left( 1 - \frac{y}{H} \right) \\
C_2 &= \frac{x}{W} \left( 1 - \frac{y}{H} \right) \\
C_3 &= \left(1 - \frac{x}{W} \right)\frac{y}{H} \\
C_4 &= \frac{x}{W} \cdot \frac{y}{H}
\end{aligned}
$$
Note: $C_1 + C_2 + C_3 + C_4 = 1 $
Then the projection on the screen is at:
$$ P = C_1 \cdot P_1 + C_2 \cdot P_2 + C_3 \cdot P_3 + C_4 \cdot P_4 $$
This is a first order approximation that should be adequate for your needs.  Especially if the TV image is a relatively small portion of the screen.
To locate the corners, a rough approach is to scan horizontally looking for black streaks.  Match them row by row.  The endpoints will define your TV edges.  Use linear regression to fit a line.  Next scan vertically in the same manner and do the same.  Calculate the intersections as your corners.  They need not be integer, and will be more accurate if you leave them fractional.
This is the easy way.  Depending upon what you are projecting will likely work fine.  The more proper way is to find the inverse mapping, then for every pixel on the screen, find the location in your source rectangle and do a weighted average there using the same splitting formula there to find the weights of the corners.


Negative values of the FFT I try to describe a simple signal in the frequency domain with the use of FFT function:

from scipy.fftpack import fftfreq, rfft

def spectrum(sig, t):
    """
    :param sig: signal.
    :param t: time scale.
    :return:
    """
    f = fftfreq(sig.size, d=t[1]-t[0])
    y = rfft(sig)
    return f, y

Signal:

...and it's FFT:

Why I'm getting negative values of the FFT? Do I even use the FFT function well? I want to check "how many" frequencies are in the shown above signal. Thank you in advance.

The rfft function returns complex values ordered in specific way. In order to retrieve the amplitude of your DFT you must take the absolute value of it. Easiest way to do that is to call y_a = np.abs(y).

I also noticed in your case that you are using fftfreq function to obtain corresponding frequency vector. In case of rfft you must use the rfftfreq. This will eliminate this nasty wrap-around on your plot. Also you should not forget about amplitude normalization by half of your window energy.

Here is the complete example:

from scipy.fftpack import rfftfreq, rfft
import numpy as np
import matplotlib.pyplot as plt

def spectrum(sig, t):
    f = rfftfreq(sig.size, d=t[1]-t[0])
    y = rfft(sig)

    return f, y

if __name__ == "__main__":
    N = 1000
    t = np.linspace(0, 4*np.pi, N)
    x = np.sin(10*np.pi*t)+np.random.randn(N)

    f, y = spectrum(x, t)

    plt.subplot(211)
    plt.plot(t,x)
    plt.grid(True)
    plt.subplot(212)
    y_a = np.abs(y)/N*2 # Scale accordingly
    plt.plot(f, y_a)
    plt.grid(True)
    plt.show()


Non-negative matrix factorization for audio separation - why does it work? Non-Negative Matrix Factorization aims to factorize a matrix $\mathbf V$ into the product of two matrices, $\mathbf V = \mathbf W\mathbf H$, where $\mathbf W$ represents a set of basis vectors and $\mathbf H$ their activations or weights. Applied to audio source separation, $\mathbf V$ is the spectrogram (which is known to us), and $\mathbf W$ is the basis non-negative spectra, $\mathbf H$ is their temporal activations.

This seems to be me to be able to express a spectrogram as a sum of different frequency bins and their magnitudes - which would help to filter or separate different pitches (seemingly doing the job of a normal band-pass filter?). I do not understand how this would actually separate, say, a piano and a singer singing the same or similar pitch.

Recall that the columns of $\mathbf{W}$ can be thought of as "basis" vectors (or elements of a dictionary - the building blocks of any signal) and elements in each column in $\mathbf{H}$ give the corresponding weights (that vary over time). This allows us to decompose the spectrogram based on not just frequency components but also temporal onset information; so it does more than what a simple bandpass filter or comb filter would do.

Non-negative matrix factorization of a spectrogram will not magically separate a piano and a singer singing in the same pitch. However it gives a useful approximation in terms of a sum of weighted basis vectors which (hopefully) can be split up into contributions from different sources because it's unlikely for each source to occupy exactly the same frequency bin at exactly the same time instant.

More concretely, let $\mathbf{V}$ be $M\times N$, $\mathbf{W}$ be $M \times K$ and $\mathbf{H}$ be $K \times N$. So we have $M$ frequency bins, $N$ time samples and $K$ decomposed components. If $\mathbf{w}_i$ are the columns of $\mathbf{W}$ and $\mathbf{h}_i$ are the rows of $\mathbf{H}$ we can write:
$$
\mathbf{V} \approx \sum_{i=1}^K \mathbf{w}_i \mathbf{h}_i^T.
$$
If we know there are only two sources in the recording (piano and singer) we can try to split up their contributions by choosing subsets of the columns of $\mathbf{W}$ and corresponding subset of the rows of $\mathbf{H}$. So if $S \subset \{1,\ldots,K\}$ then the part contributed by the piano can be denoted by:
$$
\mathbf{V}_{\mbox{piano}} = \sum_{i \in S} \mathbf{w}_i \mathbf{h}_i^T
$$
and the part contributed by the singer is:
$$
\mathbf{V}_{\mbox{singer}} = \sum_{i \in \{1,\ldots,K\}\setminus S} \mathbf{w}_i \mathbf{h}_i^T
$$

In reality, we will likely end up with a decomposition which never achieves this separation exactly. That is, there will be $\mathbf{w}_i$'s that have contributions from both the singer and the piano making it difficult to separate the two.

Here's a Python notebook showing this procedure for a mixture of drums and guitar: http://nbviewer.jupyter.org/gist/ingle/93de575aac6a4c7fe9ee5f3d5adab98f (Or if that doesn't work, here: https://gist.github.com/ingle/93de575aac6a4c7fe9ee5f3d5adab98f)

Note that the NMF algorithm only generates a decomposition. It cannot select subsets of $\{\mathbf{w}_i\}$ corresponding to each source. In the Python example, there's a manual step of figuring out which $\mathbf{w}_i$'s sound most like the guitar v/s drum. One can perhaps automate this step by noting that the drum's $\mathbf{w}_i$ vectors have more stuff in lower frequencies.

Analyzing each column (frame) of the spectrogram using a bank of bandpass filters is another way of decomposing the spectrogram. However, note that the decomposition generated by NMF is low rank i.e. parsimonious. In the Python example, it was much easier to manually select subsets of 16 columns of $\mathbf{W}$ corresponding to the two sources. With a bank of bandpass filters we would have had to turn many more knobs (# filters, locations of pass-bands for each frame) and the number of parameter combinations can grow very quickly.

References:

https://ccrma.stanford.edu/~njb/teaching/sstutorial/part2.pdf

http://musicinformationretrieval.com/nmf.html


Detect When Someone Has Stopped Talking What's an extremely efficient way of detecting when a user has stopped talking into a microphone? When I use systems like Siri it can detect when I've stopped talking almost immediately, even when there's background noise. My initial guess was to get the average volume of the background noise first but it seems there isn't much time for that as Siri starts listening straight after the button is pressed. The method doesn't have to be accurate in what is being said, just recognize that something is being said. What algorithms should I look into?

Assume noise is not a serious issue in your problem. I guess you can get pretty clear speech signal. If you have speech recognition part implemented in your system, I think you should be able to take advantage over the language model in your recognition system. According to the transition probabilities, you shall get some confidence to say at what moment someone stopped talking.

In case you donot know the language model. Here is a simple example. Suppose I want to introduce myself. I will say "Hello[verb] everyone [object]. My name [subject] is [verb] Bob [object]. I [subject] came [verb] from ... [object]". The probability for someone to stop talking right after a subject or verb is much lower than after an object. If you can recognize speech, then you shall be able to take advantage over these language models.

The other thing that is worthy to try is to see whether you can find some stop patterns when someone is talking. What I mean is that from word to word, sentence to sentence, everyone might have his/her stop patterns. If you could successfully retrieve some patterns for example duration, then they shall help you differentiate whether the current stop is alike previous stops or not.


Pulse shaping with root-raised cosine filter for a non-Gaussian channel I am here to check whether my line of thinking is correct.

Root-raised cosine (RRC) filtering is one way to minimize intersymbol interference (ISI). Hence, we apply matched RRC filters on both sides of the communication link.
We can adjust the roll-off factor $\beta \in \left(0,1\right]$ to trade-off the bandwidth (BW) for a better signal quality at the receiver. The borderline cases of $\beta$ are $\beta=0.01$ and $\beta=1$, where we "sacrifice" 0.01 BW and 0.5 BW, respectively.
Now, this is usually when the impulse response function of a channel is assumed to have a normal distribution (or similar-ish). In other words, the noise is linear (i.e., AWGN)
If this is no longer the case and the impulse response of the channel does not follow normal distribution, would RRC still be a viable solution? - I guess not but please correct me if I am wrong.

There are many related questions on this site; make sure to review them. To summarize:

Correct.
The roll-off factor $\beta \in [0,1]$. If $R_p$ is the pulse (or symbol) rate, then the pulse bandwidth is $B=(1+\beta)R_p/2$. When $\beta=0$, the pulse is a sinc and the bandwidth is the smallest possible. For $\beta>0$, the factor $1+\beta$ is the "excess" (not "sacrificial") bandwidth.
The AWGN channel does not distort the input signal, only white Gaussian noise is added to it. The density of the noise is Gaussian (normal), not the channels' impulse response. The noise, being a signal, cannot be linear or non-linear (only systems have this property).
RRC and other Nyquist pulses can avoid ISI in many different channels. For example, in  a frequency selective channel, after linear equalization the received signal has no ISI. RRC pulses remain the most popular in practice, in all kinds of channels.


correcting wobbling motion in sinogram I have measured some sinograms of a specific sample. However, there is some wobbling motion that disturbes the resulting images. Since the red upper and lower lines (indicated with black arrows) are expected to be straight without occurring wobbling motion, I want to use this information to correct my images. Is there any method that is able to straighten the lines and transform the image accordingly? I use matlab as programming language but I would also like to hear about general ideas.

The lower image is an approach to correct the image using additional information about the image (although it does not seem to work very good it might be a better starting point as the upper image)

The sinogram does make sense if the object is deformable, although if it was sagging under the force of gravity, it should be registering maximum deviation when viewed at 90 degrees (gantry 0 degees is "North"). This one seems to be registering maximum deviation at approximately 135 degrees.

If the glass tube can bend then it can compress as well. If we forget true orientation for a bit and focus on the top sinogram, then the top part of the tube is compressing while the bottom part is "expanding" to create the bend.

In the bottom image, you seem to be correcting in a "naive" way. That is by shifting columns up or down (vertically) by some amount dictated by the curvature of the walls of the tube.

However, the object is deformable. Which means that the slice that is vertical to the gantry's plane of rotation is NOT vertical to the plane of deformation. In other words, shifting up or down distorts the image because it picks up adjacent "voxels" that do not belong to the slice you are trying to correct for. Another way to look at this is that you are decreasing the resolution at the compressed areas and increasing it at the expanded areas.

This last point is important when applying corrections because it will insert some error that is proportional to the elasticity of the tube / specimen. That is, geometrically, the objects might reconstruct at more accurate locations but their density "measurements" might be slightly "blurred".

If you absolutely have to correct the sinogram then you would be looking at image warping and specifically polynomial (or spline) warping.

The approach is to fit a deformable rectangular fine grid over the deformed object and then REMAP the deformed grid on to a parallel grid using interpolation (to fill in those gaps because of curvature I am hinting above). In other words, you will be scanning the (deformed) sinogram, left-right, top-bottom ALONG the curve of deformation and each one of these points will be REMAPPED on to the straight grid locations. (Something similar to this)

This assumes that the contents of the tube are homogeneous (which is probably not true). If they are not, then their deformation depends on their elasticity, the tube deforms differently than the specimen and while you will be correcting one, you will be distorting the other.

Anyway, what you can try to do is to fit a polynomial on the top boundary (corresponding to angle) and at each point use vertical (to the slope of the polynomial) segments (corresponding to y) to recover the warping polynomials.

For more information regarding the background of polynomial image warping please see this link (Sections 2, 2.6) and this link. For more information on using MATLAB to recover and apply these transforms please see this link and more specifically this link (or possibly this link too).

Finally, an idea would be to set the tube upright and rescan but that orientation has its own set of distortions and compromises.

Hope this helps.


What is the LOW frequency resolution rule analogous to Nyquist? If I'm analyzing a signal in the frequency domain, I know about the well known Nyquist criteria that the sample frequency must be > 2x of the highest component present in the signal.

However, there are very few references related to the lowest frequency I can analyze.    This problem comes up in a lot of natural data where we have limited windows of observation because the instruments to do so were only invented or in wide spread use in the last 50-100 years, and natural cycles can be decades, centuries, or longer.

One would think that you could analyze a component where at least one period was in the sample window, but I find because of required windowing to prevent edge effects that the signal often disappears (it depends on the phase).  I then thought it was a mirror of Nyquist and I could resolve any signal with at least 2 periods in a window, but actual analysis shows that to not be true when there's a another overlying signal whose period is similar.

I found through experiment that I needed at least 5 periods to resolve a signal where there's another signal near the same frequency.

Is there a formal way of representing this problem or a definitive reference I can use?

Would the same issue apply to wavelet analysis in the same manner?

The "limit" depends on the signal-to-noise ratio, and the nearest interfering frequency(s).

If the nearest interfering frequency carrier is tuned 1% away, you might need (not just 6, but) 250 cycles or more to resolve the two separated frequency peaks, depending on the window and the S/N.  Whereas, a single pure unmodulated sinusoid in zero noise might be determined in as little as 3 or 4 non-aliased samples.

Also, there may be better low frequency estimators than depending on just a single non-rectangularly-windowed peak DFT magnitude bin, as using just that magnitude bin might be throwing some information about the phase of interference from the (nearby) conjugate negative frequency image that could be useful in improving frequency estimation and/or detection.


Spectrogram of 493 signals at a time I have been working with some seismic data and recently submitted a paper to a journal.

My seismic data has 493 traces: each trace are 8300 seconds long with sampling frequency of 2000 Hz.

One of the reviewers is asking me to generate a frequency spectrogram at the time= 2000s for all signals. My knowledge about the spectrogram is that you can only do it for one signal and there should be a time window to run the frequency analysis on. How can I do it?

Here are more detail about my problem:
These 493 traces recorded vibration around a cable in a oil producing horizontal well. The traces are 16ft away from each other along the well. The idea of the publication was to assess the spectrogram for traces  on a specific location that was the target for hydraulic fracturing. The figure below is the spectrogram of one of the traces (signals) close to the target zone. Note that there are 492 other traces and all of them have the same number of samples, frequency, and length. All these traces recorded vibration for the entire 8300secon. The question is how I could generate spectrogram of all 493 traces for example at 2000 second.
Here is the wording of the reviewer: "I welcome the addition of a spectrogram showing many channels at one particular time interval."
The paper is accepted, But I think it's better to address the reviewer's comment. I will appreciate your name in my acknowledgment.

"I welcome the addition of a spectrogram showing many channels at one
  particular time interval."

If you don't want to include it, or don't have enough space, it is always worthwhile to put them in a supplementary file, attached to the journal, or separately on the Internet. I suggest the reviewer is looking for an onset phenomenon, and it would be easy to generate a few crops in  time interval $[1500-2500]$ side-by-side, potentially for only a tenth of your data (about 50 signals).


what is the difference between the convolution and differentiation over image in image processing? In image processing the differentiation and convolution are the terms that are used interchangeably.

What is the difference between applying convolution and differentiation over the image?

How we perform the differentiation over the image(Since in order to differentiate we need to represent it as equation that is differentiable How we do it?) ?

In this context, convolution refers to the operation of LTI filtering performed on the image by the filter impulse response. Filtering can have many purposes such as blurring, sharpening, noise reduction etc.

Certain applications requires that you compute an approximation to the derivative of the image data. This can be accomplished by filering the image with a specific filter kernel (impulse response) that's some sort of a high pass characteristics.

Therefore, computing the (approximate) derivative of an image can be accomplished by LTI filtering with a highpass impulse response, which refers to a convolution operation. This is the connection in between.


Lowpass Hann filter in Python I'm trying to design a FIR lowpass filter using Hann windowing in Python but I don't know where to start from... I have fc1,fc2, ft, rp and rs as parameters. This is my Python code:

            import myDSP
        import matplotlib.pyplot as plt
        import scipy.signal as signal
        import numpy as np

        def plotInTime(x,fs):
t=np.arange(0,np.size(x,0)/fs,1/fs)
plt.figure()
plt.plot(t,x)
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.grid(True)

def plotInFrequency(x,fs):
N=int(np.size(x,0)/2)
if np.size(x,0)==1:
X=np.fft.fft(x,axis=1)
else:
X=np.fft.fft(x,axis=0)
X=np.abs(X)
X=X[:N]
f=np.arange(0,fs/2,fs/2/N)
plt.figure()
plt.plot(f,X)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Magnitude')
plt.grid(True)
        fs, x = myDSP.readWav('Ring05c.wav') #is 22050
        myDSP.plotInTime(x,fs)
        plt.title('Original Signal')
        myDSP.plotInFrequency(x,fs)
        plt.title('Spectrum for the Original Signal')

        fc1=1200
        fc2=3000
        fc= np.arange(fc1,fc2+1)
        ft=150
        rp=0.1
        rs=40
        t = np.linspace(-1, 1, 201)
        #x = (np.sin(2*np.pi*0.75*t*(1-t) + 2.1)+0.1*np.sin(2*np.pi*1.25*t + 1) +0.18*np.cos(2*np.pi*3.85*t))
        wp=(fc-ft/2)/(fs/2)
        wp = wp[0]
        ws=(fc+ft/2)/(fs/2)
        ws = ws[0]
        L, Wn = signal.cheb1ord(wp,ws, rp, rs);
        b, a = signal.cheby1(L, rp, Wn,'low')
        f,H=signal.freqz(b,a,worN=256,plot=None,fs=fs)

        y=signal.lfilter(b,a,x,axis=0)
        myDSP.plotInTime(y,fs)
        plt.title('Filtered Cheby1 Signal')

        myDSP.plotInFrequency(y,fs)
        plt.title('Spectrum for the Filtered Cheby1 Signal')
        y=y.astype('int16')

        w_fir = signal.firwin(65, fc/(fs/2), window='hann')
    [w,h] = signal.freqz(w_fir)
    plt.plot(w, 20*np.log10(abs(h)))
    plt.plot(w, np.angle(h))

You should read the documentation when using new functions: scipy.signal.firwin

You either provide pass_zero option or two cutoff frequencies. So in your case you don't have to use this option.
The first argument corresponds to the length of the filter. So it has to be an odd integer.
You either provide the sampling frequency or express the cutoff frequencies as relative frequencies (in your case something like fc1/fs).

In summary, for a a filter of order n, you can get your coefficients as

signal.firwin(n + 1, [fc1, fc2], fs=fs, window='hann')


What is "hard quantization" strategy? I am working on classification and several times I encountered with this term.
What is hard quantization strategy?
What does it differ from soft approach?

Suppose you are given a data point $D$ and want to classify it into one of $N$ possible classes, $C_i$ for $i=1,2,\ldots, N$.

Hard classification chooses one of the $C_i$ to represent the classification of $D$.

Soft classification can be done in a number of ways, but the usual (for me) is to generate a weight $w_i$ for each of the $C_i$ which indicates the probability that the data point $D$ is part of the class $C_i$. This means that
$$\sum_{i=1}^N w_i = 1$$.

I haven't heard the idea associated with quantization of one dimensional signals, but I've seen it used in relation to vector quantization which can be thought of as a simple classification technique.

This paper seems to have an example that matches what you say in your comment.

The idea is that, instead of each code word in the dictionary having the same weight, each code word has a weight dependent on some criterion or criteria: how far it is from other code words; how frequently it is likely to appear, etc. etc.

The picture below gives two examples (green and light blue) of data points that are better coded using a soft technique than a hard one.


How to prevent Octave Jumping in ACF of Speech Signal? I'm working with speech signals and my aim is to estimate the fundamental frequency
$\ F_0$   of this signal often called as "pitch".

The main idea is taking small blocks of the speech signal such that stationary can be assumed. Then calculating the autocorrelation function (ACF) of this block of a speech signal and finding index the global maximum of the ACF (except at zero) which refers to fundamental frequency.

But in the text it is stated that :

  The global maximum might not be at the lag corresponding to the true
  fundamental frequency but can possibly be an integer multiple of it.
  Due to this, the maximum can jump in consecutive frame between lags
  corresponding to multiples of T0 leading also to jumps in the
  F0-estimate. These effects are called octave-jumps.

My questions arise at this point: How does octave jumps occur? What is the possible reason? I know that ACF is a periodic function since original time sequence is periodic and in my opinion, this period equals to the block length of the original speech signal we are working on. When I investigate the interval, first period of the ACF, how can I decide whether the maximum is refers to the pitch or it is a maximum shifted from the consecutive period (block)? How can ı prevent from this effect?

Octave jumps down can occur if there is more noise (or other modulation) around a sub-harmonic of the pitch frequency, than at pitch.

Octave jumps upwards can occur if the timbre contains strong overtone or harmonic sequences and weak or nearly missing spectrum at or around the fundamental pitch frequency.  Noise, or a non-flat frequency response in the acoustic environment or channel, can also add power to an overtone peak.

Windowing effects in the ACF can add noise which can increase the probability of the above estimation errors.

Humans might be infering octave by using using preceding or following pitch, timbres, and transients to error correct spurious jumps in the time-local spectrum.


Convolution engine for FIR and IIR equal? I have a beginners question. I just installed Audiovero Acourate. This software create FIR filters for room correction of my loud speakers. The filters are .wav files which then should be used by a convolution engine, like jconvolver or brutefir on Linux.

Now, I decided to try to setup a convolution engine on Windows instead, since there are several pain points (software and driver / hardware) on Linux. I now use Equalizer APO on Windows for the convolution.

My question is: Is there any difference between convolution engines for IIR and FIR filters? And is there a difference in sound quality using either convolution engine?

I could not find much info about the convolution engine within Equalizer APO. I‘m just interested in best audio quality.

Is there any difference between convolution engines for IIR and FIR filters?
  And is there a difference in sound quality using either convolution engine?

A convolution engine performs convolution the discrete version of which is:

$$y[n] = \sum_{k=0}^{|h|-1}x[n-k] \cdot h[k]$$

Where:

$n$ denotes discrete time
$|\cdot|$ denotes length of sequence (so $|h|$ means, the length of the impulse response).
$y$ is the output
$x$ is the input
$h$ is impulse response

Notice here that this formula is giving you $y[n]$. This means that to calculate one time instance for the output you need to take into account $|h|-1$ time instances in the past of the signal ($x[n-k]$ for $k=0$ is the present value).

Both Finite Impulse Response (FIR) and Infinite Impulse Response (IIR) filters can theoretically be expressed via convolution.

But, in the case of IIR filters, $|h|=\infty$. Therefore, to produce a given $y[n]$, you would have to wait forever. What is really done in practice given an IIR filter is to truncate its impulse response when this has decayed sufficiently. The alternative (still in the case of IIR filters) is to implement them in one of the standard forms. In that case, you are implementing a recursive difference equation that produces one $y[n]$ as a combination of a finite set of past values of both its input and output. As a side note, the fact that the filter now has feedback from its output back to its input is what makes it IIR as there are now conditions by which numbers can circulate forever in that loop without necessarily decaying to zero.

Therefore:

When referring to a "...convolution engine" we are talking about FIR filters.
Because of the fact that you would have to truncate the IIR impulse response to express it via convolution, yes the result can be different.

  I could not find much info about the convolution engine within Equalizer APO. I‘m just interested in best audio quality.

Equalizer APO can do both types of filters (look into filters/) plus special implementations for some standard filters. For example, if you look into filters/IIRFilter.cpp and specifically the function IIRFilter::process(), you will see exactly the two iterations that apply coefficients to an input buffer, an output buffer and finally produce an output value. And if you look into filters/ConvolutionFilter.cpp you will notice that it performs convolution via FFT (faster) and libHybridConv (see libHybridConv/libHybridConv.c, hcProcessSingle()).

In conclusion: "Best audio quality..." depends on the workload of your machine. A sufficiently long IIR filter would not sound too different than an FIR one. And at the same time, long filters take more time to calculate. This is further worsened if you want to apply the filter in more than one channels (and I am not talking multiple speakers here, just going stereo means applying one convolution per channel. You should also take into account the comments about the sound quality being subjective here :) ).

However, you might be able to take the results of your room's impulse response, prioritise which "problems" you want to correct and then use a simple graphic equaliser to correct those problems in particular.

You can run a simple graphic EQ, whether in Equalizer APO or any other VST host and then re-route your audio source to be going first through the VST filter and then to the output.

Hope this helps.


What are the proper pre-processing steps to perform Independent Component Analysis? What are the proper steps for preprocessing my waveforms in order to perform an independent component analysis (ICA) later? I understand the how, though further explanation of that doesn't hurt, but I'm more interested in the why.

Independent component analysis (ICA) is used to separate a linear mixture of statistically independent and most importantly, non-Gaussian† components into its constituents. The standard model for a noise-free ICA is

$$\mathbf{x}=\mathbf{As}$$

where $\mathbf{x}$ is the observation or data vector, $\mathbf{s}$ is a source signal/original components (non-Gaussian) and $\mathbf{A}$ is a transformation vector that defines the linear mixing of the constituent signals. Typically, $\mathbf{A}$ and $\mathbf{s}$ are unknown.

Pre-processing

There are two main pre-processing strategies in ICA, namely centering and whitening/sphering. The primary reasons for pre-processing are:

Simplification of algorithms
Reduction of dimensionality of the problem
Reduction of number of parameters to be estimated.
Highlighting features of the data set not readily explained by the mean and covariance.

From the introduction of G. Li and J. Zhang, "Sphering and its properties", The Indian Journal of Statistics, Vol. 60, Series A, Part I, pp. 119-133, 1998:

  Outliers, clusters or other kind
  of groups, and concentrations near curves or non-ﬂat surfaces are probably the
  important features that interest data analysts. They are, in general, not obtainable through mere knowledge of the sample mean and covariance matrix. In
  these circumstances, it is desirable to separate oﬀ the information contained in
  the mean and the covariance matrices and forces us to examine aspects of our
  data sets other than those well-understood natures. Centering and sphering is a simple and intuitive approach that eliminates the mean-covariance information
  and helps to highlight structures beyond linear correlation and elliptic shapes,
  and therefore is often performed before exploring displays or analyses of data
  sets

1. Centering:

Centering is a very simple operation and simply refers to subtracting the mean $\mathbb{E}\{\mathbf{x}\}$. In practice, you use the sample mean and create a new vector $\mathbf{x}_c=\mathbf{x}-\overline{\mathbf{x}}$, where $\overline{\mathbf{x}}$ is the mean of the data. Geometrically, subtracting the mean is equivalent to translating the center of coordinates to the origin. The mean can always be re-added to the result the end (this is possible because matrix multiplication is distributive).

2. Whitening:

Whitening is a transformation that converts the data such that it has an identity covariance matrix, i.e., $\mathbb{E}\{\mathbf{x}_c\mathbf{x}_c^T\}=\mathbf{I}$. Normally, you work with the sample covariance matrix,

$$\widehat{\mathbf{\Sigma}}=C.\mathbf{x}_c\mathbf{x}_c^T$$

where $C$ is just my lazy placeholder for the appropriate normalization factor (depending on the dimensions of $\mathbf{x}$). A new whitened vector is created as

$$\mathbf{x}_w=\widehat{\mathbf{\Sigma}}^{-1/2}\mathbf{x}_c$$

which will have a covariance of $\mathbf{I}$. Geometrically, whitening is a scaling transformation. Here is a small example in Mathematica:

s = RandomReal[{-1, 1}, {2000, 2}];
A = {{2, 3}, {4, 2}};
x = s.A;
whiteningMatrix = Inverse@CholeskyDecomposition[Transpose@x.x/Length@x];
y = x.whiteningMatrix;
FullGraphics@GraphicsRow[
  ListPlot[#, AspectRatio -> 1, Frame -> True] & /@ {s, x, y}]

The first plot is the joint density of two uniformly distributed random vectors, or the components $\mathbf{s}$. The second shows the effect of multiplying by a transformation vector $\mathbf{A}$. The square gets skewed and scaled into a rhombus. By multiplying with the whitening matrix, the joint density is back to a square which is slightly rotated from the original.

Because of the whitening transformation, in the new system that is being solved, i.e. $\mathbf{x}_w=\mathbf{A}_w\mathbf{s}_w$, $\mathbf{A}_w$ is an orthogonal matrix. This can be easily shown:

$$
\begin{align}
\mathbb{E}\{\mathbf{x}_w\mathbf{x}_w^T\}&=\mathbb{E}\{\mathbf{A}_w\mathbf{s}_w(\mathbf{A}_w\mathbf{s}_w)^T\}\\
&=\mathbf{A}_w\mathbb{E}\{\mathbf{s}_w\mathbf{s}_w^T\}\mathbf{A}_w^T\\
&=\mathbf{A}_w\mathbf{A}_w^T=\mathbf{I}
\end{align}
$$

where the last step follows because of the statistical independence of $\mathbf{s}_i$ The orthogonality condition means that there are only about half as many parameters that need to be estimated. (Note: Although this is true in this case and in my example, $\mathbf{A}$ need not be square to begin with).

If, after the transformation, there are eigenvalues close to zero, then these can be safely discarded as they are just noise and will only hamper the estimation due to "overlearning".

3. Other pre-processing

There might be other pre-processing steps involved in certain specific applications that are impossible to cover in an answer. For example, I've seen a few articles which use the log of the time-series and a few others that filter the time-series. While it might be suited for their particular application/conditions, the results don't carry over to all fields.

†I believe it is possible to use ICA if at most one of the components is Gaussian, although I can't find a reference for this right now.

Why is it called "sphering"?

This is probably well known, but just as a fun fact, sphering comes from the change in the structure of covariance matrices in the case of Gaussian components from an $n$-dimensional  hyper ellipsoid to an $n$-dimensional sphere due to whitening. Here's an example (use the same code as above, but replace {-1,1} with NormalDistribution[])

The first is the joint density for two uncorrelated Gaussians, the second under transformation and the third is after whitening. In practice only steps 2 and 3 are visible.


Fourier Transform for frequency analysis I applied Fourier transform on a sound (PCM generated digitally with sin function) and frequency detection is accurate. However, if I play the sound and record it, and then apply a Fourier transform to detect frequencies, my results are nowhere close to original results. The window size I have chosen for Fourier transform is 30 samples. ( Sampling frequency is 44100 Hz and 30 samples represent a symbol. Each symbol frequency is shifted by 1470 Hz.). Is my approach right ?

That's to be expected. If the input of your FF transform is 30 (real) points, you also get 30 real amplitudes out, for 30 different frequency bins. Each bin is (44100/2)/30 Hz wide, or about 700 Hz. That is quite wide. Add to that the effect of spectral leakage, which may quite well extend for 5 bins in either direction. 10 out of 30 bins will contain the ground frequency. Any higher harmonics will add to that.

Normally, people use much longer FFT's. 256 is quite a common choice, and 1024 isn't unheard of either.


Downsampling fft results to a graphics display I'm new to dsp.  I've got a complex signal that was sampled at 192kHz then an fft of 16384 bins created.  This now needs to be plotted on a graph for a real time spectrum display.  The problem is if you want to plot the whole spectrum consisting of a 192kHz span, you will be plotting 16384 points on a display that many times will be 1600 pixels wide for example.  The result is you get an ugly mess on the screen.  So far I have been downsampling the signal by just plotting every 8th point. The problem is though, when I look at other spectrum displays that also perform 16384 fft's, they show much more data.  I have effectively dropped some on the peaks that the other software doesn't.
So the question is, how does one do the downsampling properly?
Thanks

Ok, here is the answer...hopefully this will help someone else.  Read this thesis:

http://skemman.is/stream/get/1946/15343/37285/3/SS_MSthesis.pdf

This gives a very easy and comprehensive overview of how to do just what I needed. I do believe this particular case is different than the usual downsampling. It fits the situation perfectly.
The thesis describes several methods of downsampling for display purposes. For example the practice of just dropping every nth sample, performing an average over the bucket of samples corresponding to the downsample ratio and more advanced methods.  The one I will be looking at is called Largest Triangle Three Bucket.
Basically you divide your original samples into buckets consisting of samples equal to your downsampling factor.  For example if you have a sample set of 10,000 points and you want to only display 2000, you will divide your original samples into 2000 buckets containing 5 sample each.
You then rank every point in the bucket by calculating the area of a triangle it forms
with the selected point in the last bucket and the average point in the next
bucket.
This seems to work really well in preserving the visual representation of data.


highpass filter design in Matlab I hope to design 1st order highpass filter from transfer function. In the example of 1st order lowpass filter, I first get the coefficients of numerator and denominator in the variable 'b' and 'a'. In case of first order filer,

b = 1/tau;
a = [1, 1/tau];
h = freqs(b, a, w) % LPF
[digital_b, digial_a] = bilinear(b, a, fs) % Analog to Digital transformation.

The order is like above. First, we obtain the coefficients of the transfer function and then make LPF by the 'freq' function in MATLAB. After that, I can convert the analog filter to a digital filter by using the 'bilinear' function.

However, in case of HPF, I cannot transform it using the 'bilinear' function. When I run the 'bilinear' function in MATLAB for 1st order HPF, it gives the following error

Numerator cannot be higher order than denominator.

I used the code below to make digital HPF.

b = [1, 1/tau];
a = 1;
h = freqs(b, a, w) % LPF
[digital_b, digial_a] = bilinear(b, a, fs) % Analog to Digital transformation.

What is the problem? and how should I design HPF on MATLAB?

Transfer function of a 1st order HP filter is:

$$H(s)=\frac s{s+wc}=\frac{{\displaystyle\frac s{wc}}}{1+{\displaystyle\frac s{wc}}}\;,where\;wc=2\mathrm\pi\ast\mathrm{fc}$$

&

a = [1 w0];
b = [w0 1e-48]; % ! two coefficients


frequency domain zero padding in OFDM I have a matlab code for OFDM system, in this code the auther added a zero padded subcarriers before the IFFT stage, in the middle of the raws, and added zero padded guard interval after IFFT, my question is that why i need to add zero padded subcarriers and what is the relation between the number of null subcarriers and sampling frequency??

Zero padded OFDM (i.e.) adding zeros after IFFT is used to combat Inter symbol interference, suppose that there is multi path in the system that spans over L samples then atleast L zero samples are prefixed at the transmitter after IFFT  so that these can be thrown away at the receiver and then recover ISI free IFFT frame / OFDM symbol.

The advantage of zero padding also has to do with power saving at the transmitter, the transmitter is simply switched off during the transmission of zero padded L samples duration. This has the advantage of saving power when compared to cyclic prefix OFDM which should then also transmit the cyclic prefix.

Zero padded OFDM though can leave to noise boosting and this is a disadvantage when compared to cyclic prefix OFDM. The equalization also would be non linear in case noise boosting is to be avoided. There is a lot of literature on zero padded OFDM you could refer to that as well.

Regarding the null sub carriers at the ends, these are implemented to reduce the ACLR and make the design of anti aliasing filters easier ( the cut off can be less sharp) hence save cost


NLMS algorithm for a MISO structure I am trying to implement an NLMS algorithm for a multi-input single-output(MISO) structure.

We take a reference signal x, then we made a new set of P input signals from it as follows: x_op (k) = x(k)^p. k denotes the k-th sample of our reference signal x.

For the case where P = 1, our adaptive filter is just a vector of length N. But for P > 1, we have a matrix H (N x P) containing all the adaptive filters. (where N is the length of each filter).

Now I am confused as what would be dimenstion of the estimated output at each iteration (for each sample). For the case where P = 1, we can simply write:

y_hat = h_hat' * xk; where xk has the last N samples of x.

and y_hat is then a scalar. What about this case when H is a matric and not a vector anymore?

The way I understand it, y_hat for this case is no longer a scalar. But if it's not a scalar, then how I am supposed to define the error for each sample k, in order to write the adaption law equation?

I can provide you coding example, but it solves more generic problem than you posted. It is for general dynamic system that converts two inputs x(t), y(t) into single output z(t) and defined by Urysohn integral equation. Hammerstein is particular case of Urysohn and linear system is particular case of Hammerstein. The code and some explanation is here:
http://ezcodesample.com/NAF/MISO.html
http://ezcodesample.com/NAF/Wiener.html
http://ezcodesample.com/NAF/reallife2.html

And most generic case when dynamic system is an operator that converts fragment of input x(t) from t - T to t into scalar z(t) http://ezcodesample.com/NAF/kolmogorov.html


Calculating the SNR from the power spectrum. Is my suggestion relevant? I have a signal that is sampled at $20 \mathrm{Hz}$. It has almost stationary behavior, with a main frequency around $1\ \mathrm{Hz}-1.5 \ \mathrm{Hz}$. This frequency range is to be expected and is of interest for further processing.  Unfortunately, I don't know the noise present at the signal in order to calculate the SNR of the signal as a way of assessing its quality.

I was thinking that since I know a priori the frequency range of the fundamental frequency of interest that should exist in the signal, I could calculate the SNR based on the power spectral density (PSD) of the signal. Moreover, after PSD calculation using Welch's method, I could integrate within a frequency range that corresponds to my main frequency (assuming that within this frequency range there is one main lobe) and divide that to the integration of the rest of the calculated spectrum as an SNR metric:

$$\mathrm{SNR}_\textrm{metric}=\frac{\int_{n_{f_1}}^{n_{f_2}}p_{xx}(n)*f(n)dn}{\int_{n_\textrm{f=0}}^{n_{f=f_1}}p_{xx}(n)*f(n)dn + \int_{n_{f=f2}}^{n_{f=\textrm{flast}}}p_{xx}(n)*f(n)dn}$$

where $p_{xx}$ the PSD, $f(n)$ the corresponding frequencies vector based on Welch's PSD calculation method, $n$ is the index of the frequencies and $f_1$, $f_2$ is the bandwidth at the PSD that surrounds the fundamental peak of the frequency of interest that should be around 1.0 to 1.5 Hz when the lobe attenuates completely.

Would you consider this approach valid and why?
If you have advice for a different method to calculate the SNR let me know!

I would like to thank everyone in advance.

Assuming the noise is evenly distributed among all frequencies, this approach makes sense in principle. However, some corrections:

You should not integrate over $pxx(n)f(n)dn$, since n is already the frequency. I think, you mean $\Delta_f(n)$, which is the frequency step between two bins. Assuming $\Delta_f(n)$ is equal and independent of n, then you can even cancel it in the equation.
You need to divide the noise and signal components by the respective bandwidth of the integral. Otherwise you will calculate more noise, when the bandwidth increases.
currently, you calculate (signal+noise)/noise. Though, you want to calculate signal/noise, i.e. (signal+noise-noise)/noise


fourier series fitting matlab I am using cftool in Matlab to fit time series values to Fourier model with 8 terms:

 General model Fourier8:
 f(x) =
           a0 + a1*cos(x*w) + b1*sin(x*w) +
           a2*cos(2*x*w) + b2*sin(2*x*w) + a3*cos(3*x*w) + b3*sin(3*x*w) +
           a4*cos(4*x*w) + b4*sin(4*x*w) + a5*cos(5*x*w) + b5*sin(5*x*w) +
           a6*cos(6*x*w) + b6*sin(6*x*w) + a7*cos(7*x*w) + b7*sin(7*x*w) +
           a8*cos(8*x*w) + b8*sin(8*x*w)

But I don't understand how Matlab calculate a and b coefficients. I found in literature that I can use formulas:

ak=2/n * sum ( x(t) cos (2*pikt/n))

bk=2/n * sum ( x(t) sin (2*pikt/n))
 where n is number of terms in time series, k is index of current term k=1,...n, sum is with t,

but I can't get values which Matlab showed me for these coefficients.

Is Matlab General model Fourier8- formula for trigonometric Fourier series? Because in literature the first term in formula is a0/2, here is a0? Also Matlab value for a0 is not mean value. How to calculate this correct?

The cftool is going to use a curve fitting procedure.  Some time of non-linear regression.  In a general sense, it is a trial and error method:
1) guess a value for the coefficients 2) calculate the function with the coefficient, and determine the error between the calculated function and the original data, 3) try a new value for the coefficients, that is slightly different than the original guess. 4) go back to step 2, and keep going until the error is sufficiently small.

cftool is a very general tool that can be used for all kinds of problems, not just fourier series. cftool has no idea that you are trying fit a fourier series. In a general sense, it really is no different than if you just tried a whole bunch of combinations until you found one that worked, it's just that because it is a computer it can try many thousands of combinations in a second.

The formulas you showed, on the other hand, are the actual expression for calculating a fourier series.  Depending on your data, cftool will probably never match these exactly, although it might get close.  The reason it doesn't match is that there may be many different combinations that match the data fairly well, and cftool cannot tell the difference between them.  One way to help cftool get closer would be to give it more data to work with (i.e. a longer sample).

Depending on what your goal is, you might just want to use the closed form expression for a fourier series.  that will be the most accurate, and probably quicker to compute anyway.


Given two waves, 441 Hz and 700 Hz, how many samples long is each wave's period if the sample-rate is 48 kHz? The title says it all.

I have two waves, Master and Slave.  Master is at 441 Hz and Slave at 700 Hz.  The sample-rate is 48 kHz.

I want to know how to calculate the period of each wave in terms of samples.  How many samples long is each wave's period?

Please show me how to calculate this

Thanks

Neither of these tones is an exact number of samples at 48 kHz, but you can calculate the periods as follows:

$\frac{48000}{441} = 108.84$ samples

$\frac{48000}{700} = 68.57$ samples


Inverse Laplace Transform A system given by $\frac{s-1}{(s+1)(s-2)}$ has to be inverse transformed so that it is anticausal and nonstable. The answer given is $h(t)=-\frac{1}{3}(2e^{-t}+e^{2t})u(-t)$

Why the minus sign at the beginning?

First you have to remember a Laplace transform property:

$$ e^{a t} u(t) \longleftrightarrow \frac{1}{s-a} ~~~,~~~ \mathcal{Re}\{s\} > \mathcal{Re}\{a\}   \tag{1} $$

$$ -e^{a t} u(-t) \longleftrightarrow \frac{1}{s-a} ~~~,~~~ \mathcal{Re}\{s\} < \mathcal{Re}\{a\}   \tag{2} $$

This property states that a given Laplace transform $H(s)$ can correspond to multiple inverse transforms, depending on the region of convergence ROC.

Therefore given $H(s) = \frac{1}{s-a}$ , you can find two possible inverses as $x(t) = e^{at} u(t)$ or $x(t) = -e^{at} u(-t)$ depending on whether the ROC is to the left or right of the pole location. Note that one of them is causal and the complementary is anti-causal.

Now, given your transfer function

$$ H(s) = \frac{s-1}{(s+1)(s-2)} = \frac{2/3}{s+1} + \frac{1/3}{s-2} ,  \tag{3} $$

it has two poles at $s = -1$ and $s=2$. There will be three possible ROC's with three different inverses :

$$
\begin{align}
\mathcal{Re}\{s\} < \mathcal{Re}\{-1\} & \implies h(t) = -\frac{2}{3} e^{-t}u(-t) - \frac{1}{3} e^{2t} u(-t)  \tag{4} \\
\mathcal{Re}\{-1\} < \mathcal{Re}\{s\} < \mathcal{Re}\{2\}  &\implies h(t) = \frac{2}{3} e^{-t}u(t) - \frac{1}{3} e^{2t} u(-t)  \tag{5} \\
\mathcal{Re}\{s\} > \mathcal{Re}\{2\} & \implies h(t) = \frac{2}{3} e^{-t}u(t) + \frac{1}{3} e^{2t} u(t)  \tag{6}
\end{align}
$$

The impulse response in (4) is anti-causal and unstable.

The impulse response in (5) is non-causal and stable.

The impulse response in (6) is causal and unstable.


Confusion related to PN sequence : Terminology and application in statistical signal processing A Pseudo-random Noise (PN) sequence is a sequence of binary numbers, e.g. ±1, which appears to be random; but is in fact perfectly deterministic. The sequence appears to be random in the sense that the binary values and groups or runs of the same binary value occur in the sequence in the same proportion they would if the sequence were being generated based on a fair "coin tossing" experiment. In the experiment, each head could result in one binary value and a tail the other value.

Consider a model, $$ x(t) = h^T u(t) +w(t)  \tag{1}$$
where $u$ is the input and $w$ is Additive White Gaussian noise.

In Book : Fundamentals of Statistical Signal Processing by Steven Kay, Chapter 4 : Eq(4.21) finds the variance of the estimated coefficients when the input is PRN as

$$var(\hat{h_i} )= \sigma_w^2/Nr_{uu}[0] \tag{2}$$ under the assumption of the input to be a PN .
where $\sigma_w^2$ is the variance of the measurement noise; $\sigma_u^2$ is the variance of the input
The MVU estimator using a PRN sequence is mentioned as $\hat{\theta }= (H^TH)^{-1} H^T x$

It is stated that for the deterministic sequence $u$ the $H^TH$ becomes a symmetric Toeplitz autocorrelation matrix. For $H^TH$ to be diagonal we should use a PRN sequence as input which yileds $H^TH = Nr_{uu}[o]I$.
These expressions can be found in document Pg 7 Section 4.1.3

My problem is that I am confused if PRN means binary or $\pm1$ sequence or only $iid$ real valued White Gaussian noise. If PRN is different from PN or binary, then what will be the estimates and the CRLB?

If the output $x[\cdot]$ of a discrete-time linear system is related to the input via
$$x[n] = \sum_{i=0}^m h[i]u[n-i] \tag{1}$$
where $u[\cdot]$ is a periodic sequence with period $N \gg m$, then $x[\cdot]$ is also a sequence with period $N$ and
the
periodic cross-correlation function $R_{x,u}[\cdot]$ between $x[\cdot]$
and $u[\cdot]$ is the periodic convolution between the sequence $h[\cdot]$
and the periodic autocorrelation function $R_{u,u}[\cdot]$. That is,
$$R_{x,u} = x \star u = (h * u) \star u = h * (u \star u) = h* R_{u,u}
\tag{2}$$
  If
$u[\cdot]$ is a binary pseudonoise (PN) a.k.a. pseudorandom noise (PRN) sequence that takes on values $\pm 1$, then
$$R_{u,u}[n] = \begin{cases}N, & n \equiv 0 \bmod N,\\
-1, &\text{otherwise}\end{cases}\tag{3}$$
which is essentially a periodic digital impulse function,
 and by ignoring that $-1$ as inconsequential compared to $N$, we get from $(2)$ that
$$N h[n] = R_{x,u} \implies h[n] = N^{-1} R_{x,u}[n].\tag{4}$$
If we include additive white Gaussian noise $w[n]$
of variance $\sigma^2$ on the right side of $(1)$,
then $(2)$ becomes
$$R_{x+w,u} = (x + w)\star u = (h * u + w) \star u = h * (u \star u)
+ w \star u= h* R_{u,u} + \hat{w}
\tag{5}$$
where $\hat{w} \sim \mathcal N(0, N\sigma^2)$. Thus, the measured
cross-correlation function is now a Gaussian random variable
with mean $Nh[n]$ and variance $N\sigma^2$.  The estimator for
$h[n]$ is thus just $N^{-1}$ times the measured value of the
cross-correlation function, and its variance is thus
$\displaystyle \frac{\sigma^2}{N}$, This matches up perfectly with
Kay's result
$$\operatorname{var}(\hat{h}[n]) = \frac{\sigma_{w}^2}{Nr_{u,u}[0]}$$
since we have assumed that $u[n] \in \{+1, -1\}$ and hence
$r_{u,u}[0] = 1$: all the $u[n]$, regarded as random instead of
pseudorandom, have variance $1$.


How to uniquely identify phase difference between harmonics? Let's say I have a multi-frequencies signal of 440Hz and 880Hz, slightly out of phase.

I can take a chunk of the signal and calculate phase difference of the two frequencies with FFT bins. However, since frequencies are different, the phase difference will differ depending upon where I took the chunks, despite they are taken from the same signal.

I wonder whether there is a way to uniquely identify phase difference between harmonics, with arbitrary chunks? (so that I can identify harmonic signals.)

i am assuming this is a periodic (or quasi-periodic) function.  also seems like the signal is a musical tone.  periodic signals have a period, $P$, and the reciprocal of that period is the fundamental frequency; $f_0 \triangleq \frac{1}{P}$

or, as angular frequency; $\omega_0 = 2 \pi f_0 = \frac{2 \pi}{P}$ .

the periodic signal is $x(t)$ where

$$ x(t+P) = x(t)  \qquad \forall t \in \mathbb{R} $$

it can be represented as a Fourier series

$$ \begin{align}
x(t) \ &= \ a_0 \ + \ 2 \sum_{k=1}^{\infty} a_k \cos(k \omega_0 t) - b_k \sin(k \omega_0 t) \\
 &= \ a_0 \ + \ 2 \sum_{k=1}^{\infty} a_k \frac{e^{j k \omega_0 t} + e^{-j k \omega_0 t}}{2} \ - \ b_k \frac{e^{j k \omega_0 t} - e^{-j k \omega_0 t}}{2j} \\
 &= \ a_0 \ + \ \sum_{k=1}^{\infty} (a_k + jb_k) e^{j k \omega_0 t} \ + \ (a_k - jb_k) e^{-j k \omega_0 t} \\
 &= \ \sum_{k=-\infty}^{\infty} c_k \ e^{j k \omega_0 t} \\
\end{align} $$

where

$$ c_k =
\begin{cases}
a_{-k} - jb_{-k}, & \quad \text{for } k < 0 \\
a_0, & \quad\quad k = 0 \\
a_k + jb_k, & \quad\quad k > 0
\end{cases} $$

and, going in the other direction,

$$ \begin{array}{lcl}
a_0 & = & c_0 \\
a_k & = & \Re\{c_k \} \quad\quad \text{for } k>0\\
b_k & = & \Im\{c_k \} \quad\quad\quad k>0
\end{array} $$

note that for real $a_k$ and real $b_k$, then

$$ c_{-k} = c_k^* = \text{"complex conjugate of } c_k \text{ "}$$

it's easier to get $c_k$ from $x(t)$ and then compute $a_k$ and $b_k$ from the $c_k$ which is:

$$ c_k = \tfrac{1}{P} \int\limits_{t_0}^{t_0+P} x(t) \, e^{-j k \omega_0 t} \, \mathrm{d}t \qquad \text{for any } t_0 \in \mathbb{R}$$

you can uniquely define the phase of overtones (or harmonics apart from the fundamental) with reference to the fundamental if there is any energy at the fundamental.  (it is not always the case that a periodic signal has non-zero energy in its own fundamental frequency.)

now it's also true that:

$$ \begin{align}
x(t) \ &= \ \sum_{k=-\infty}^{\infty} c_k \ e^{j k \omega_0 t} \\
 &= \ a_0 \ + \ 2 \sum_{k=1}^{\infty} a_k \cos(k \omega_0 t) - b_k \sin(k \omega_0 t) \\
 &= \ a_0 \ + \sum_{k=1}^{\infty} r_k \cos(k \omega_0 t + \phi_k) \\
\end{align} $$

where

$$ \begin{align}
r_k &= 2 \big| a_k + j\,b_k \big|  \\
 &= 2 \sqrt{a_k^2 + b_k^2}  \\
\\
\phi_k &= \arg \{a_k + j\,b_k \}  \\
 &= \operatorname{atan2}(b_k, a_k)  \\
\end{align} $$

now, the phases of all of the harmonics are relative to some time $t=0$, but that reference time can be defined in such a way that the phase of the fundamental is zero.

let

$$\begin{align}
 \tilde{x}(t) \ &\triangleq \ x(t - \tfrac{\phi_1}{\omega_0}) \\
 &= \ a_0 \ +  \sum_{k=1}^{\infty} r_k \cos \big(k \omega_0 (t - \tfrac{\phi_1}{\omega_0}) + \phi_k \big) \\
 &= \ a_0 \ +  \sum_{k=1}^{\infty} r_k \cos(k \omega_0 t + \phi_k - k\phi_1) \\
 &= \ a_0 \ +  \sum_{k=1}^{\infty} r_k \cos(k \omega_0 t + \tilde{\phi}_k) \\
\end{align} $$

where $ \tilde{\phi}_k \triangleq \phi_k - k\phi_1 $.

you can see that $\tilde{\phi}_1 = 0$, that we defined our origin of $t$ so that the phase of the fundamental (when $k=1$) is 0.  then all of the other harmonics have their phases referenced to that "first harmonic" or fundamental.  the only problem is that the amplitude of that fundamental must not be zero.  $r_1 \ne 0$ so that the reference phase, $\phi_1 = \arg \{a_1 + j\,b_1 \}$, can be determined.


Writing a Discrete Fourier Transform program I would like to write a DFT program using FFT. This is actually used for very large matrix-vector multiplication (10^8 * 10^8), which is simplified to a vector-to-vector convolution, and further reduced to a Discrete Fourier Transform.

May I ask whether DFT is accurate? Because the matrix has all discrete binary elements, and the multiplication process would not tolerate any non-zero error probability in result. However from the things I currently learnt about DFT it seems to be an approximation algorithm?

Also, may I ask roughly how long would the code be? i.e. would this be something I could start from scratch and compose in C++ in perhaps one or two hundred lines? Cause actually this is for a paper...and all I need is that the complexity analyis is $\mathcal{O}(n \log n)$ and the coefficient in front of it doesn't really matter :) So the simplest implementation would be best. (Although I did see some packages like kissfft and FFTW, but they are very lengthy and probably an overkill for my purpose...)

When applied correctly, the DFT yields exactly the same result as convolution, it is no approximation (neglecting rounding errors due to limited precision of floating point calculations, of course).

To calculate the circular convolution of vectors $x$ and $y$ (both of length $N$):
$$
x \circledast y = \operatorname{IDFT}_N\left[\operatorname{DFT}_N[x]\cdot\operatorname{DFT}_N[y]\right]
$$
To calculate the linear convolution first append $N-1$ zeros to both vectors yielding the zero-padded vectors $\tilde x$ and $\tilde y$. Then:
$$
x * y = \operatorname{IDFT}_{2N-1}\left[\operatorname{DFT}_{2N-1}[\tilde x]\cdot\operatorname{DFT}_{2N-1}[\tilde y]\right]
$$

To implement the [I]DFT I highly recommend the FFTW library. I promise that it will take you less time to learn its usage than to write your own implementation. For the special case of $x$ and $y$ containing only $1$ and $0$ their might be a more efficient algorithm, though. You would have to look into the Cooley-Tukey algorithm that is used for the FFT and see if you can simplify it for your needs. Whether it gets simpler will also depend on how you define computational complexity - if you define it in terms of additions and multiplications a custom FFT implementation probably won't buy you much. But if you go down to the hardware implementation level you can save some multipliers in the first butterfly stage by replacing them with multiplexers. Anyway, it is well known that the FFT algorithm complexity is of order $N\log N$ as you say. I would like to think that you don't have to show this again.

If you'd like to play around with IFFT and FFT I recommend Octave (an open source Matlab clone).


Choosing Codes or sequences with excellent Auto-Correlation properties The Auto-Correlation function of Walsh-Hadamard codewords does not have a good characteristics.

It can have more than one peak and thus, the Walsh-Hadamard codes do not have the best spreading behavior or correlation property.

The cross-correlation function of the Walsh-Hadamard codewords can also be non-zero for a number of time shifts and unsynchronized users can interfere with each other.

My Questions:

Which other codes or sequences show :

A good auto-correlation attribute?
Good Cross-correlation Property ?

There is a vast literature on this subject. In particular, deleting one bit
(the leading bit in one specific implementation) from each of the Walsh-Hadamard sequences and applying a permutation to the remaining $2^n-1$ bits will result
in $2^n$ sequences of length $2^n-1$ that consist of

(i) the all-zeroes sequence

(ii) the $2^n-1$ cyclic shifts of a pseudonoise (PN), or maximal-length linear  feedback shift register (LFSR), sequence

A PN sequence has ideal autocorrelation properties, but the cross-correlation
between a sequence and its cyclic shift has a peak value that is the same as
the autocorrelation peak.

For more than you probably want to know, I refer you to the paper
D.V. Sarwate and M.B. Pursley, "Cross-correlation properties of pseudorandom and related sequences," Proc. IEEE, vol.68, pp.593-619, May 1980.  It includes a
detailed discussion of the Gold sequences mentioned by JasonR as well as the
small sets and large sets of Kasami sequences that you will encounter in some
applications.


Design $ {L}_{2} $ Norm Optimal Infinite Impulse Response (IIR) Filters It is widely known that matching a FIR filter of fixed length to a band model is an unconstrained QP-problem. The MATLAB function firls() implements a solution to this problem. Basically, one symbolically evaluates integrals of the (Z-domain) transfer function over the unit circle, getting astonishingly simple expressions for the match quality. Additionally, that method is stable and will not fail, other than the Parks McClellan Filter Design Algorithm.

Is it possible to fit IIR filters to such a model in the $ {L}_{2} $ Norm Sense?
Several problems arise:

While the cost function (fit quality?) can be written down, is it not a simple quadratic anymore, but a fraction of two quadractics.
You appear to have to optimize nominator and denominator at the same time.
It might be hard to constrain the denominator to "minimum phase", so that the resulting filter is stable.

I am not looking for gradient descent or other local optimization methods, as the cost function might have several local optima. An idea that crossed my mind would be that it appears possible to pose the FIR problem so that it yields the optimal coefficients for an IIR all-pass.

As you have already pointed out in your question, it is not possible (without using optimization methods) to compute an exact L2 solution for the frequency domain design problem of IIR filters due to the non-linear relationship between the filter coefficients and the error function. There is, however, a method which can come close and which transforms the problem to a linear one: the Equation Error Method. Instead of defining the error measure as
$$E_0=\sum_k\left |H_k -\frac{B_k}{A_k}\right|^2\tag{1}$$
(where $k$ is the frequency index, $H_k$ is the complex desired frequency response at frequency $\omega_k$, and $B_k$ and $A_k$ are the numerator and denominator polynomials, respectively, also evaluated at the frequency point with index $k$), one defines an error measure
$$E_1=\sum_k\left |H_kA_k -B_k\right|^2\tag{2}$$
Minimizing $E_1$ is a linear problem in the filter coefficients. You get an overdetermined system of linear equations which can be solved (in the $l_2$ sense) by solving a set of linear equations. (Of course, the number of frequency points must be greater than the number of filter coefficients). The result of minimizing (2) is identical to solving a weighted $l_2$ problem with weight function $|A_k|^2$:
$$E_1=\sum_kW_k\left|H_k -\frac{B_k}{A_k}\right|^2\quad\textrm{with}\quad W_k=|A_k|^2$$
There are two (related) problems with this method:

The specification must not only include the desired magnitude response but also the desired phase. If the phase specification is chosen in a way that does not fit the chosen filter order and the general properties of IIR filters, the approximation error will be large, and the filter might be unstable (which brings us to the next point).

Stability is not considered in the design process. Depending on the specification, the best approximation could either be stable or unstable.

Of course, poles outside the unit circle can be reflected inside the unit circle without affecting the designed magnitude response, but the approximation error might still be large, because an unstable filter indicates that the specification is not well suited to the chosen filter order. So it takes quite some experience to choose an appropriate desired phase response.
An in-depth explanation of the equation error method can be found in the book Digital Filter Design by Parks and Burrus. You can find a good overview in Julius O. Smith III - Introduction to Digital Filters - With Audio Applications - Equation Error Formulation.


Visualizing aliasing in Matlab I am trying to verify for myself that aliasing actually makes it impossible to distinguish the "real" and the aliased frequency. What I have done it set the sampling rate to 20 Hz and created two sine wave with 2 Hz and 18 Hz. As far as I understand the 18 Hz signal should alias to a 2 Hz signal when the sampling rate is 20 Hz.

Matlab code:

t = [0:0.05:1];      //20Hz sampling
a = sin(2*pi*2*t);   //2Hz sine wave
b = sin(2*pi*18*t);  //18Hz sine wave
plot(t, a);
hold on;
plot(t, b, 'red');

The plot looks like this:

I can see that both sine waves are equal, but 180 degrees out of phase. I wonder if there is something to it or if it's just Matlab drawing through the points in opposite ways?

I think you should plot something like:

t = [0:0.05:1];      %20Hz sampling
a = sin(2*pi*2*t);   %2Hz sine wave
b = sin(2*pi*18*t);  %18Hz sine wave
plot(t, a, 'bo');
hold on;
plot(t, b, 'ro');

T = [0:0.001:1];     %1000Hz sampling frequency
A = sin(2*pi*2*T);
B = sin(2*pi*18*T);
% plot for 1000Hz sampling frequency
plot(T, A, 'b');
plot(T, B, 'r');

legend({'2Hz @ 20 Hz f_s', '18Hz @ 20 Hz f_s', '2Hz @ 1000 Hz f_s', '18Hz @ 1000 Hz f_s'})

That allows to understand which exact samples of your sinusoids you are taken at a given sampling frequency:


Understanding Arrangement of second order sections When splitting higher order recursive filters into more convenient 2nd order sections (biquads), in matlab, we usually use zp2sos to convert. In document, there's a couple of lines , which I am quoting below,

[sos,g] = zp2sos(z,p,k,''order') specifies the order of the rows in sos, where 'order' is
'down', to order the sections so the first row of sos contains the poles closest to the unit circle
'up', to order the sections so the first row of sos contains the poles farthest from the unit circle (default)

Here i can't remember that I had read something relevant to arranging the biquads on the basis of poles location.
Now Can anyone explain me,
why should we arrange the biquads in a specific manner while cascading them ?
What that arrangement is to deal with the poles being closer or farthe from unit circle, as Matlab documents says?
What specifics need to be considered while ordering the biquads to get higher order filters.
It'll be so nice to illustrate by using some case study, for example, designing the A weighting filter (8th order IIR), and converting it ino biquads, and arranging them in a specific manner which is under discussion.
Thanks in advance.

Different orderings of second-order blocks have different behavior with respect to quantization noise and overflow. Usually it is best to avoid very peaky transfer functions in second order blocks, because these peaks may cause overflow and quantization noise amplification. This is why it is common to pair poles and zeros that are close to each other, so the zeros cancel part of the peaks caused by the poles. The pole-zero pairing determines the second-order sections, but not their order. Concerning the order, there is always a best choice minimizing the quantization noise, but this order needs to be determined individually for the filter under consideration. However, there are two types of ordering that work well in most cases. Depending on the way scaling is performed, either ordering the sections with increasing or decreasing pole radii is best. In Matlab, these two common ways of ordering are referred to as 'up' and 'down', respectively. Have a look at this document for more detailed information.


Confusion in N of N-point-FFT On what basis, N of N-point-fft is chosen? What can go wrong if I don't select it to be  power of 2? How is this N related to sampling frequency?

Firstly, you need to know that the FFT computes the DFT (discrete Fourier transform) in an efficient manner. Hence, the output of an N-point FFT and N-point DFT are exactly the same.

Now, especially, if N is a power-of-two, the FFT can be calculated very efficiently. Though, if N can be factorized into rel. small prime numbers (2, 3, 5), still the calculation can be done very efficiently. Strictly speaking, the Radix-2 FFT only works, if N is power-of-two. However, in common implementations in software, N can be any number and the implementation chooses the appropriate algorithms for most efficient calculation.

Hence, if N is power of two, the calculation by FFT is just a bit faster, compared to when N would be a large prime number.

Regarding sampling frequency: If $F_s$ is your sampling frequency and $T$ is the duration of the FFT window, then $N=TF_s$. It's hard to answer more detailed what you mean by relation to sampling frequency, because I dont know exactly what you are asking.

Regarding your comment: The DFT maps a subset of the frequency range $[0,\dots,F_s]$  to the $N$ DFT-bins, which are indexed by $k=0,\ldots,N-1$. Since all frequency bins in the DFT output have the same distance to each other (i.e. they are linearly spaced), you immediately end up with $f=\frac{k}{N}F_s$, i.e. the $k$th bin corresponds to the physical frequency $f=\frac{k}{N}F_s$.


Separate speech signals that are bandlimitted I have 2 kind of signals (refer figures) in a folder. One is band limited to 4khz and another is band limited to 8khz, both are with sampling frequency 16khz. I need to separate them. Any idea??

As mentioned in comments that the 4Khz and 8 Khz bandlimited files are kept separately but in same folder and problem is to segregate them into separate folders.

One can use scripts written in Matlab or Python for doing this, procedure could be as follows :

Read each file (this could be easily automated in matlab/python) from folder.
split each file into frames, take FFT of frames.
for the 4 Khz band limited files the magnitude of upper half of FFT coefficients will be very less, but for 8Khz magnitude of upper half of FFT coefficients will be significantly high.
the above logic could be automated by keeping fixed threshold on magnitudes of FFT coefficients.
based on above threshold, files could be moved to separate folders.


Image region-based representation I have a homogeneity image that includes two regions $\Omega_1$ and $\Omega_2$. How to represent the image by decomposition of two regions?

For example, the famous way is representing by the mean feature. Image $I$ can be  represented as follows:

$$I=u_1c_1+u_2c_2$$

$c_1$ is mean of region $\Omega_1$

$c_2$ is mean of region $\Omega_2$

$ u_1 =
  \begin{cases}
   1 & \text{if is $\Omega_1$ }\\
   0       & \text{otherwise}
  \end{cases}
$

$ u_2 =
  \begin{cases}
   1 & \text{if is $\Omega_2$ }\\
   0       & \text{otherwise}
  \end{cases}
$

Do you know another way to represent image $I$ without using mean feature. Thank you so much

As the question, I will show a simple example as following

Your example seems to imply you are talking about mutually exclusive, non-overlapping regions.

In this case, given N such regions, use an image where each pixel has enough values (depth), to represent N regions. Just set the pixel value to the index of the region in that position.

You could think of the result image as a sum of N binary images, each multiplied by the $i$ subscript of the corresponding $\Omega_i$.

Alternatively, if you are asking about how to segment such an image into 2 (or more) regions, the you should use a connected component algorithm, which is common is most image processing libraries.


opencv otsu binarization not working I tried to apply opencv OTSU binarization on the following image, but it generate a black image.

Here is the following python code:

import cv2

image = cv2.imread('test_otsu.png', 0)
ret, thresh = cv2.threshold(image, 0,250,cv2.THRESH_BINARY+cv2.THRESH_OTSU)
cv2.imwrite('temp.png', thresh)

Did I miss anything in my code? Please help me how to get correct binarization of the given image.

the problems seems to be with the white border of the image you are using. remove it manually and test again.


Could modulation/coding schemes be generated from an AI? Warning: I'm a bit of a noob with respect to this entire field, but I took some classes on it (a few years ago) and found it fascinating.

Anyhow, as far as I understand it, modulation is the process of converting data to a signal, using quite a fair bit of creativity in the process. For example QPSK can be presented in the complex domain in the following way:

Timing diagram for π/4-QPSK. The binary data stream is shown beneath the time axis. The two signal components with their bit assignments are shown the top and the total, combined signal at the bottom. Note that successive symbols are taken alternately from the two constellations, starting with the "blue" one. (from wikipedia)

The resulting signal is also very confusing:

Timing diagram for π/4-QPSK. The binary data stream is shown beneath the time axis. The two signal components with their bit assignments are shown the top and the total, combined signal at the bottom. Note that successive symbols are taken alternately from the two constellations, starting with the "blue" one. (from wikipedia)

The actual question

Could an AI be used to generate better modulation/coding schemes? Perhaps kind of like the following AI "learned" to walk: https://www.youtube.com/watch?v=gn4nRCC9TwQ, we could create a "game" (really an advanced physics simulation) in which the AI has to modulate and demodulate a signal and try to achieve a very high data rate while keeping bit errors as low as possible in different conditions. These new standards would probably be even harder to explain to humans but more efficient.

Any help is welcome:

Has this already been done?
Has there been any research on this? If yes, how much more efficient could generated processes be?
Is this not possible yet because it would be complicated to create transmitters and receivers for the process?

As per my knowledge, one of the recent papers on this topic can be found here, where the authors used machine learning algorithms to generate the optimal constellation shaping. Surprisingly, the results were in line with an old work: F. R. Kschischang and S. Pasupathy, "Optimal nonuniform signaling for Gaussian channels," in IEEE Transactions on Information Theory, vol. 39, no. 3, pp. 913-929, May 1993.

Moreover, autoencoders, neural networks, reinforcement learning and deep learning are being applied for end-to-end learning (especially in the case of wireless communications), where the mathematical model of the wireless communication channel is not very accurate. One such article can be found here.

Also, for the case of polar codes, deep-learning-based decoders are very popular nowadays at least in academia.


Coordinate system in digital image processing I am trying to understand image coordinate system where I read that it is reversed.
So if $f(x,y)$ is actually $f(y,x)$ where $f(\cdot, \cdot)$ is the image. I am a bit confused and would appreciate help in clarifying certain points.

In digital image processing, an image is a matrix of rows ($m$) and columns ($n$) denoted by $[m,n]$. When showing as an image output would the rows become the $Y$ axis and and column the $X$ axis? Just wanted to confirm this as sometime it is confusing. For example, consider a $3 \times 5$ matrix:

[ 30   40  20  30  0 60
  10   20  15  10  0 20
  20   5   10  30  0 10 ]

where the elements represent the pixel intensities. If I want to access the element 60, how would I specify? Would row, $m = 1$ and column, $n = 5?$ Is row the $X$ or $Y$ axis?

In Maths, if I were to plot in a graph the coordinate $(5,1)$ then 5 refers to the $x=5$ (plotted on horizontal $X$ axis) and $y=1$ (plotted on vertical $Y$ axis). However, in image representation the coordinate $(5,1)$ means go to row 5 and column 1 ie., vertical axis is the row, $m=5$ and the horizontal is the column, $n=1$? Is there a reversal going on or am I misunderstanding something. What is the convention?

So in general, in image Columns --> X axis and Rows ---> Y Axis? But in graphs it is the opposite?
UPDATE:

Consider a matrix
A = [1,2,3
     4,5,6
     7,8,9]

On the graph the coordinate point is $(3,2)$. Here $x=3$, $y=2$. If I want to access element at $A(3,2)$ then this means row (horizontal) = 3, column (vertical) = 2 and the element at $A(3,2) = 2$.
In images does the vertical axis become the rows and the horizontal the columns or is it the same as that of matrices, except for the top to bottom traversal constraint? Something seems to be changing, that is the confusing part.

Different operating systems, standards, and graphical softwares implement distinct (but quite related) conventions for representing, storing, manipulating, or displaying two-dimensional data on computers. The distinction is about the orientation of the axes.
Basically we have 4 main object types :

$f(x,y)$ , x-y continuous, is a function of mathematical relevance, such as an analog image.
$f[n,m]$, n-m integers, is a sequence obtained by sampling $f(x,y)$, or  synth-generated.
$F(i,j)~$, i-j integers, is a matrix to represent the data in its rows-i and columns-j.
$f[i,j]~~$ , i-j integers, is an array to store the data used in a computer program.

The arguments $x,y,n,m,i,j$ can also be replaced with $x_1,x_2,n_1,n_2,i_1,i_2$, the subscripts indicating the dimension order of the variables, such as the first or the second dimension.
A sequence $f[n_1,n_2]$ can represent a 2D digital image, or a mathematical function defined over the coordinate system of axes $n_1,n_2$. Its samples are placed in the correct elements of the associated matrix $f(i,j)$ for the operation to produce the desired effect at the output.
MATLAB fundamentally uses MATRIX based dimension ordering. The first sample A(1,1) is at the top-left corner, first dimenison $i$ is along vertical-down the rows, and the second dimension, $j$, is along the columns (horizontal-right).
This is equivalent to placing a coordinate system of $n_1 ~, ~n_2$ with its origin (0,0) at the element A(1,1); the first axis $n_1$ points along the first dimension $i$ (vertical down); and the second axis $n_2$ points along the second dimension $j$ (horizontal right) of the matrix A(i,j).
This mapping is obtained when you rotate the conventional coordinate system ($n_1$ horizontal-right, and $n_2$ vertical-up, origin at bottom-left) by 90 degrees clockwise with respect to its origin (0,0) and placing the origin at the element A(1,1) of the matrix.
In a parallel fashion, the 2D-DFT (discrete Fourier transform) $F[k_1,k_2]$ of the sequence $f[n_1,n_2]$ is computed by the statement:  F = fft2( f ) which places the first dimension, $k_1$, along the rows $i$, and the second dimension, $k_2$, along the columns $j$ of the output matrix $F(i,j)$ to represent the DFT sequence $F[k_1,k_2]$. Which is also alligned with the first and second dimenions of the input matrix $f(i,j)$ that represents the seqeunce $f[n_1,n_2]$.
Below is an oldskool discussion of a few mapping modes, and functions in MATLAB related with data orientation. As long as mappings are used consistently, they all yield the same results when interpreted correctly.
% SEQUENCES, BMP IMAGES, MATRICES and MATLAB FUNCTIONS :
% ------------------------------------------------------
%
% All data processing and display in Matlab is done via MATRICES A(i,j).
% But, theory of image processing is based on SEQUENCES f[n1,n2], F[k1,k2].
% When processing images, orientation of axes become relevant across
% functions such as CONV2(), FFT2(), IMSHOW(), STEM3(), SURF().
% And a MAPPING convention from f[n1,n2] into A(i,j) should be used.
%
%
% A matrix A(i,j) is indexed by vert rows i, and horz columns j.
% A sequence f[n1,n2] is indexed by horz-right n1, and vert-up n2.
% Which is the most typical, and natural, orientations for n1 and n2.
%
%
% The mapping convention depends on the functions being called :
%
%           1- PROCESS functions: conv2(), fft2(), filter2()
%           2- DISPLAY functions: imshow(), surf(), stem3()
%
%
% We consider following mapping modes between the samples of a sequence
% f[n1,n2] and elements of a matrix A(i,j) as follows:
%
%
% MM-0 : 90d CW ROTATED MAPPING :
% ----------------------------------------------------------------------
%  f[0,0] --> A(1,1), and "n1" grows DOWN from the top-row of A.
%
%    n1 = i-1   , n2 = j-1     ===>          A(i,j) = f[i-1, j-1]
%    i  = n1+1  , j  = n2+1    ===>        f[n1,n2] = A(n1+1,n2+1)
%
%  ---o---1--------2-----> j (n2)
%   1 | f[0,0]  f[1,0]
%   2 | f[1,0]  f[1,1]    Ex mapping of f[n1,n2] into 3x2 matrix A(i,j)
%   3 | f[2,0]  f[2,1]
%     |
%   i v
%  n1
%
%
%
% MM-1 : FLIP-DOWN MAPPING :
% ----------------------------------------------------------------------
% f[0,0] ---> A(1,1) and "n2" grows DOWN from the top-row of A.
%
%    n1 = j-1   , n2 = i-1     ===>          A(i,j) = f[j-1, i-1]
%    i  = n2+1  , j  = n1+1    ===>        f[n1,n2] = A(n2+1,n1+1)
%
%  ---o---1--------2-----3---> j (n1)
%   1 | f[0,0]  f[1,0]  f[2,0]
%   2 | f[0,1]  f[1,1]  f[2,1] Ex mapping of f[n1,n2] into 3x2 matrix A(i,j)
%     |
%   i v
%  n2
%
%
% MM-2 : BMP MAPPING :
% ----------------------------------------------------------------------
% f[0,0] is stored in A(N2,1) and "n2" grows UP from the bottom-row of A
%
%    n1 = j-1   , n2 = N2-i     ===>          A(i,j) = f[j-1, N2-i]
%    i  = N2-n2 , j  = n1+1     ===>        f[n1,n2] = A(N2-n2,n1+1)
%
%  n2 ^
%   2 | A(1,1)  A(1,2)         An example mapping into 3x2 matrix A
%   1 | A(2,1)  A(2,2)
%   0 | A(3,1)  A(3,2)
%  ---o---0--------1-----> n1 (j)
%
%
%
% Axis orientations of functions CONV(),FILTER(),FFT2():
% ----------------------------------------------------------------------
% They assume 90d ROTATED / or FLIP DOWN index mapping by default.
%
%      A(1,1)---------> n2,k2,j (horizontal - column variable)
%           |
%           |  A(i,j)
%           |
%           v n1,k1,i (vertical - row variable)
%
%
% NOTE: Strictly speaking, MATLAB does not care which mapping was used to
% generate the matrix A(i,j). Rather it treats "i" as the first dimension
% and "j" as the second dimension, and outputs acordingly.
%
%
% Axis orientation of DISPLAY functions STEM3(),  SURF()            :
% -------------------------------------------------------------------
% It produces the plots according to NATURAL X-Y orientation in which
% n1 point horizontal right, and n2 point verticcal up and (0,0) being
% at the bottom. But this requires that the MATRIX A(i,j) was filled in
% according to FLIP-DOWN mapping mode.
%
%                          A(N2,N1)
%                          f22
%                       f12   f21
%     A(N2,1) (n2-i) f02  f11  f20 (n1-j) A(1,N1)
%                      f01  f10
%            n2-axis      f00       n1-axis
%                        A(1,1)
%
%
% Image display function IMSOW(A(i,j)) assumes a BMP mapping:
% -------------------------------- --------------------------
% IMSHOW displays the matrix contents in its row-column order.
% Therefore if you want to display a sequence f[n1,n2], or F[k1,k2] using
% the IMSHOW() function, then in order to have the correct orientation
% according to natural x-y coordinates, you should use BMP mapping mode
% on the matrix which is to be displayed.
%

SUMMARY
If IMSHOW() iwill be used to display images, or their FFT results, then correct alignment with $n_1,n_2$ and $k_1,k_2$ requires BMP based mapping to be used to fill in the associated matrix. If STEM3() or SURF() will be used to get a 3D display of the sequences, or FFT results, then flipp-down mapping produces correct orientation. Other functions CONV2(), FFT2() etc., work equally well with either Rotated or Flip-Down mapping modes.


Problem to implement a sinusoid oscillator in matlab I'm pretty new in dsp and I'm having tough times trying to get the output for a simple digital oscillator. I'm starting from the recursive formula $y[n]=x[n]+a \times y[n-1]-b \times y[n-2]$, and need to figure out the values of $a$ and $b$ for the output to be completely sinusoidal at a determined frequency $\omega _0$ when an unit impulse is applied as input. I can't get the theory too much and the code I wrote just doesn't work at all. I would appreciate a lot with your help with this, cheers.

The code I've tried so far:

n=0:2*pi/200:2*pi;

w0=pi/16;

f=w0/(2*pi);

beta=1;

alpha=2*cos(w0);

num=1;

den=[1 -alpha beta];

x=[0 0 1 zeros(size(n))];

y=zeros(size(n));

y(1)=0;

y(2)=0;

y(3)=sin(w0);

for i=3:size(n)

    y(i)=x(i)+alpha*y(i-1)-beta*y(i-2);

end

y;

stem(n,y)

Thanks very much

OK, since the original code was changed a bit, here is the complete solution (with some useful modifications):

n=1000;
w0=pi/16;
beta=1;
alpha=2*cos(w0);
y=zeros(1,n);
y(1)=0;
y(2)=sin(w0);

for i=3:n
    y(i)=alpha*y(i-1)-beta*y(i-2);
end

stem(y)
% End of code


Purpose of FIR Filter Structure I read Digital Signal Processing of author G.Proakis and see this FIR Filter Structure. I understand all about this but I wonder that what is purpose of classification to 4 structures?
Here is this section:

These are just different ways to implement FIR filters. In theory these structures are all equivalent in the sense that they compute exactly the same output signal for a given input signal (under ideal circumstances, i.e., infinite precision), but in practice they behave differently. The differences are mainly the following:

their behavior when implemented with finite word-length
computational complexity
memory requirements

I'm sure that if you continue reading the book, these structures will be explained in more detail, and their relative merits will be discussed.


Is it correct to subtract two signals acquired over different times and trials to remove the common signal in them? I have two signals that represent the response of a neuron under two different conditions.

Signal 1 (S1): response to Stimulus A
Signal 2 (S2): response to Stimulus A+B

The response to stimulus A is common in both signals but I am interested in the response to stimulus B only. Because the signals are noisy, they have been acquired multiple times. S1 is acquired 60 times and S2 is acquired only 4 times.

Will averaging S1 over these 60 acquisitions and S2 over 4 (Figure 1) and then subtracting the averaged signals (S2_avg - S1_avg) give me a better representation of response to Stimulus B or should both the signals be averaged across equal acquisitions (Figure 2)?

Averaging S1 60 times really smooths the response to Stimulus A but compared to this, response to Stimulus A is not smooth in S2 as it is only averaged 4 times. So in this way, what I understand is  by averaging S1 many times might actually be a disadvantage.

  Additional info following answers: my underlying assumption is that
  the neuron behaves like an LTI system. My concern is that averaging S1
  many times smoothes it as compared to response to stimulus A in S2. As
  such, subtracting them would not really remove the common signal in
  them.

Will averaging S1 over these 60 acquisitions and S2 over 4 (Figure 1) and then subtracting the averaged signals (S2_avg - S1_avg) give me a better representation of response to Stimulus B or should both the signals be averaged across equal acquisitions (Figure 2).

Ideally, both signals should be averaged over the same amount of repetitions and the ideal here would be the biggest of the two (60). This is for each individual point estimate to "suffer" from approximately the same (and as low as possible) variance.

Potentially, you might gain some more insight if you were to use bootstraping, where you would be forming point estimates of each time instant of your time series and then examine these smaller subsets (potentially through a boxplot) to try and "guess" what is the most likely time series to be formed, give these two datasets. But the fact that averaging is suggested in the first place anyway, means that an assumption is already made for the time instance differences to be coming from a gaussian distribution.


phenomenon when window size for STFT is too short Hi I was researching on STFT when I encountered this slide online

I don't understand the points, Can anyone explain them? Thanks in advance.

It might mean that if you take a window of not enough length so as to contain at least one period of a sinusoid, and instead you are only getting a lot less inside the window of analysis (Like 10% of the period or something like that), then the fourier transform will show you the spectrum of a signal as if that little piece you got inside the window, extended periodically was the entire signal.

To know exactly what is going on in your graph, we need to know the scale of the time axis, and which signal was being analyzed, if it was a sinusoid, what was the frequency. but that might be what the book is trying to say.


Convolution in frequency domain Simple math question. The convolution theorem states that multiplication in time domain is equal to convolution in frequency domain and vice versa. There is a condition that the signal has to be properly zero padded as to not cause aliasing.

This question concerns convolution in the frequency domain. The difficulty arises from the fact that we are dealing with a complex signal which has a positive and a negative side. Let's consider two frequency domain signals, with negative frequencies included, presented in polar coordinates []:

a=[0,1(3π/4),0,1(π/2),0)]
b=[1,1,1,1,1]

So we have a sine wave and a zero-phase dirac delta. It's apparent that the time domain product is zero. But coming to this conclusion via using convolution in frequency domain doesn't seem to be so simple. I get (rectangular):

[-j,-j,-j,-j, 0,j,j,j,j,j]
=> [-2j,-2j,0,2j,2j]

Which is not correct. Even for a cosine I get a row of twos, which is not correct either (should be row of ones). Where am I going wrong?

If you have two DFTs $A[k]$ and $B[k]$ (note the correct representation of a sinusoid at DFT bin number $1$)

A = [0,-j,0,0,j]; B = [1,1,1,1,1];

with the corresponding time-domain sequences $a[n]$ and $b[n]$

a = ifft(A);    % [0, 0.38042, 0.23511, -0.23511, -0.38042];
b = ifft(B);    % [1,0,0,0,0];

then the multiplication of the time-domain sequences $c[n]=a[n]b[n]$ corresponds to the cyclic (or circular) convolution of the DFTs $A[k]$ and $B[k]$:

$$\text{DFT}\{c[n]\}=C[k]=\frac{1}{N}\sum_{n=0}^{N-1}A[n]B[k-n]_{\text{mod} N},\quad k=0,1,\ldots,N-1\tag{1}$$

where the indices are taken modulo $N$ ($N=5$ in this case), such that they remain in the range $[0,N-1]$:

$$\begin{align}C[0]&=\left(A[0]B[0]+A[1]B[N-1]+\ldots+A[N-1]B[1]\right)\cdot\frac{1}{N}\\
C[1]&=\left(A[0]B[1]+A[1]B[0]+\ldots+A[N-1]B[2]\right)\cdot\frac{1}{N}\\&\vdots\\
C[N-1]&=\left(A[0]B[N-1]+A[1]B[N-2]+\ldots+A[N-1]B[0]\right)\cdot\frac{1}{N}
\end{align}$$

It's easy to see that the cyclic convolution of $A[k]$ and $B[k]$ as given above results in the zero vector, because due to $B[k]=1$ for all values of $k$, each element of the result is the sum of all elements of $A[k]$, which is always zero.