How to test digital filters? First the question(s):

How should I write unit tests for a digital filter (band-pass/band-stop) in software? What should I be testing? Is there any sort of canonical test suite for filtering?
How to select test inputs, generate expected outputs, and define "conformance" in a way that I can say the actual output conforms to expected output?

Now the context:
The application I am developing (electromyographic signal acquisition and analysis) needs to use digital filtering, mostly band-pass and band-stop filtering (C#/.Net in Visual Studio).
The previous version of our application has these filters implemented with some legacy code we could use, but we are not sure how mathematically correct it is, since we don't have unit-tests for them.
Besides that we are also evaluating Mathnet.Filtering, but their unit test suite doesn't include subclasses of OnlineFilter yet.
We are not sure how to evaluate one filtering library over the other, and the closest we got is to filter some sine waves to eyeball the differences between them. That is not a good approach regarding unit tests either, which is something we would like to automate (instead of running scripts and evaluating the results elsewhere, even visually).
I imagine a good test suite should test something like?

Linearity and Time-Invariance: how should I write an automated test (with a boolean, "pass or fail" assertion) for that?
Impulse response: feeding an impulse response to the filter, taking its output, and checking if it "conforms to expected", and in that case:

How would I define expected response?
How would I define conformance?

Amplitude response of sinusoidal input;
Amplitude response of step / constant-offset input;
Frequency Response (including Half-Power, Cut-off, Slope, etc.)

I could not be considered an expert in programming or DSP (far from it!) and that's exactly why I am cautious about filters that "seem" to work well. It has been common for us to have clients questioning our filtering algorithms (because they need to publish research where data was captured with our systems), and I would like to have formal proof that the filters are working as expected.

Some thoughts on it at least.

First, if you can shown linearity and time-invariance and you know that the filter has the correct impulse response you are home. Given that the filter is stable. So in this case there is no need to run different other input signals and the frequency response is given based on the impulse response.

Checking impulse response is obviously quite simple. Linearity and time-invariance I guess is actually really complicated (any specific input signal you test may just work out of coincidence, but others may not). However, as the filter code in most cases shouldn't be that complicated, it should be possible to eye-ball the code and see that it does what is expected (no non-linear operators, no if-cases based on signal values etc).

Remains stability. For IIR filters this is a really complicated issue and I will not go into that discussion. For FIR filters, they will always be stable (unless they are implemented using some recursive algorithm). However, you may run into numerical issues relating to overflow (most likely not underflow) and round-off noise. You can find formal proofs that an FIR filter will not overflow (maximum output is the sum of the absolute values of the impulse response coefficients times the maximum input magnitude). For round-off errors you can find statistical expressions, but they are not of much use in formally proving anything (but rather as an argument that the expected output with "infinite precision" and the actual output still conforms).

Also not that all impulse responses etc should be evaluated using the actual (possibly rounded) coefficients not the (possibly higher precision) ones obtained from the filter design tool. This is because the filter will always realize the transfer function given by the coefficients that it use, no matter how they were derived. This may still not be the impulse response you get at the output though, as there may be round-off issues. An easy way to see that is to use a scaled impulse (i.e. not one) and note that typically not only the magnitude of the frequency-response will change, but also the shape. The filter will still implement the same transfer function, but the round-off noise added to the two different input signals are different.

Hope this gives you a bit more insight.


Uniqueness of Fourier Series Representation and the Fourier Transform of Periodic Signals If we are given a signal of the form $$x(t) = \sum_{k = -\infty}^{+\infty} a_k e^{j k \omega_0 t},$$ can we call it a Fourier Series representation of $x(t)$ right away?

Suppose we are given the signal $x(t) = e^{j2\Omega t} + e^{j4\Omega t}$, this signal can be expressed as $$x(t) = \sum_{k = -\infty}^{+\infty} a_k e^{jk\Omega t},$$ however, in this case, $2\Omega$ is the fundamental frequency of the signal, not $\Omega$ as the equation above might suggest.

I ask this because in the second edition of Signals and Systems by Alan Oppenheim, he derives the Fourier transform of a periodic signal by considering the impulse train $$X(j\omega) = \sum_{k = -\infty}^{+\infty} 2\pi a_k\delta (\omega - k \omega_0)$$
and applying the inverse Fourier transform to obtain $$x(t) = \sum_{k = -\infty}^{+\infty} a_k e^{j k \omega_0 t},$$ which he says is the Fourier series representation of the signal without further discussion. But how do we know that in this case $\omega_0$ is the fundamental frequency of the signal? Couldn't it be the case that it is a signal of the form given above, namely, $x(t) = e^{j2\omega_0 t} + e^{j4\omega_0 t}$, in which case the fundamental frequenct is $2\omega_0$? Thank you in advance.

The Fourier series, if it exists, is unique, apart from the fact that we can always use a multiple of the fundamental period $T$ as the period, or, equivalently, a fraction of the fundamental frequency $\omega_0$. The only consequence of this is that we introduce Fourier coefficients with value zero in between the non-zero coefficients. The values of the non-zero coefficients remain unchanged:

$$x(t)=\sum_{k=-\infty}^{\infty}a_ke^{jk\omega_0}=\sum_{k=-\infty}^{\infty}\tilde{a}_ke^{jk\omega_0/N},\quad N\in\mathbb{Z}^+$$

$$\tilde{a}_k=\begin{cases}a_{k/N},\quad k=lN,\quad l\in\mathbb{Z}\\0,\quad\text{otherwise}\end{cases}$$


The exact definition of dominant frequency? What is the exact definition of dominant frequency? I have googled and cannot find it.

For example, does this figure shows "two" dominant frequencies? Is dominant frequency a single maximum value? or do we call it dominant frequency because it has a distinguishable peak observed by human?

Edit: another example,

As dominant frequency it is usually meant the one that carries more energy w.r.t. all the other frequencies in the considered spectrum.

Since in the first example  the two frequencies have the same absolute value I suppose they are both dominant frequencies. However the concept is not really important in these case, since you are just looking at the combination of two pure sinusoids of equal amplitude.

In the second case, the lowest frequency is both the fundamental and the DOMINANT frequency since it has the highest peak of the considered spectrum.

Reference: https://arxiv.org/pdf/1306.0103.pdf


Discrete signals notation underlined value In a book i have there's an exercise that want to find the output of a system for a given  signal

system: $ 2y(n)=2x(n)+x(n+1)+x(n-1) $

input : $ Χ=[\underline1, -1, 2, 0, -2, 0, 1 ]$

$Χ1=Χ[0] $

What represent's the underlined value?

An arrow or an underline is used to indicate the value at $n=0$. Example,

$$y(n) = \begin{array}{ccccc}[&9 &8 &\underline7 &6 &5 &]\end{array}$$ means that,
$$\begin{array}{ccccc}&y(-2) = 9, &y(-1) =8, &y(0)=7, &y(1) = 6 &\mathrm{and} &y(2) =5\end{array}$$

In your case, it says that $$\begin{array}{ccccc}&X(0) = 1, &X(1) =-1, &X(2)=2, &\ldots\end{array}$$

Note: No underline or arrow is equivalent to underlining the first value.


harmonic product spectrum on guitar string signals i am trying to create a program that gets the fundamental frequency of the signal coming from a guitar string to identify what note is being played, is harmonic product spectrum the algorithm the most suitable in getting the fundamental frequency? is there any algorithm to use that can be used real-time?

also i am using matlab simulink to implement this algorithm to arduino

all answers are deeply appreciated, thank you in advance :)

HPS (harmonic product spectrum), ASDF, AMDF, and autocorrelation are all similar periodicity estimation methods (they use different weightings and algorithms), which might be useful for guitar note pitch estimation.  These methods have different estimation accuracy/error statistics, which may vary depending on the types of guitars, strings, microphone setup, and etc., tested.  Experiment.

HPS has one interesting advantage.  The overtone series of the lower guitar strings can be slightly inharmonic, partially due to the physics of actual guitar strings, which have non-zero thickness and stiffness.  However, there is some experimental evidence that the stretch between the inharmonic overtone series might actually correspond more closely to the typical human psychoacoustic perception of pitch, than the actual fundamental vibration mode frequency of the guitar string.

That said, naive use of these methods may be too heavyweight for an Arduino AVR processor.  Perhaps an Arduino using an ARM M0 would be more suitable.  Also, to reduce the computational load of doing a full FFT or full autocorrelation-equivalent, you may want to use a modified zero-crossing pitch estimator.  Just make sure to use a large enough set of zero crossing span pairs (not just adjacent crossings) so that you can use something like a single AMDF to pick which set of span pairs represent overtone, harmonic or noise lags, and which span pairs might represent fundamental pitch lags.  Once you validate your set of span pairs or lags, you may want to average, or median filter.


Generate Minimum distance codes I am trying to generate a set of linear block codewords of length n and have minimum hamming distance dmin. I understand Hamming codes but they are not flexible with n. One possible way I can think of is brute force search. But when trying to code this in MATLAB I am not sure how to proceed.

Google search helped me find this link where it does exactly what I want.
Code Generation Tool

Also I would appreciate any other methods/algorithms which can achieve this in a more efficient way.

You don't really put any constraints on the family of codes (not even the alphabet). For real applications, you do need to put more constraints on what you want for a code. And distance may not be so important (e.g. the minimum distance of many LDPC codes is quite low, but you need to decode well beyond distance/2 to come close to capacity on many channels).

Brute force searching a code is a bad idea, since most codes are hard to decode (you need some structure to exploit to decode efficiently). Also, computing the minimum distance is NP hard in general (Vardy, 1997).

Here is an example of a code, which has worked well for many applications (CD's, RAID arrays, TV, Voyager mission, etc.). Reed-Solomon Codes are a family of maximum distance separable (MDS) codes, i.e. for a given length n and number of information symbols k, they have the largest distance possible (n-k+1).

They are also easily decodable, say by the Peterson-Gorenstein-Zierler algorithm or Berklenkamp-Massey algorithm (more efficient). You can apply these to generalized reed solomon codes as well (See R.M. Roth, Introduction to coding theory, for details). Many other codes are related to this family. For example, BCH codes are alternant codes from a RS code.


Compare two wave files I have two .wav audio(stereo) files to be compared. The files may vary in frequency or amplitude or phase shift. How to compute the difference between the wave files in Scilab?  The expected output is that if frequency and amplitude differ between the audio files, the output must return that they are dissimilar.   Phase changes need not be compared.

Short Time Fourier Transform can be used for comparing the frequency and amplitude differences between two audio files. it can be implemented using basic windowing of audio followed by DFT operation over the windowed audio.


What is the input dimension for LMS update algorithm I just realized this forum existed, posted my question on the wrong forum.
Sorry about the formatting, please don't hesitate to ask if any part of this post is unclear. This is my first post.
I'm trying to implement a simulation of an ANC system with python, using this model here, but that's not the main point. My simulation keeps diverging, so I'm trying to figure out if it's my system is unstable or did I implement the NLMS algorithm wrong. Below is my code for the NLMS update. I'm using a source from Mathworks here.
x is input vector--------------------------------------------dimension(N * 1)
error is (desired vector - NLMS prediction)---------dimension(N * 1)
w_prev is the previous iteration's weight vector---dimension(1 * m)
def lms_update(x,error,w_prev,mu):

    ep = 0

    w_prev = np.array(w_prev)
    error = np.array(error)

    x_normal = np.dot(np.conjugate(x),np.linalg.inv((ep+np.dot(np.transpose(np.conjugate(x)),x))))

    w_update = w_prev + 2*mu*np.dot(error,x_normal)
    w_update = w_update.tolist()

    return w_update

then from above dimensions, after applying the update formula,
w(n)=αw(n−1)+f(u(n),e(n),μ)  and  f(u(n),e(n),μ)=μe(n)u'(n)
the update would be a 1 * 1 constant, and the weight would all be updated with the same value, which doesn't make too much sense to me.
So I adjust the input vector x into a (N * m) matrix, using vectors from current and previous iterations, so I can get the update to be a (1 * m) matrix. But I'm not sure if this is the correct way. From what I gathered from papers and google searches, nothing pops up indicating dimensions of the input. But I did implement a gradient descent algorithm before which input is a matrix, not a vector.

The output of an FIR filter is
$$
y(n) = \sum_{i=0}^{N-1}w_i(n)x(n-i)
$$
Say the weight vector $\mathbf{w}(n)$ has a length of $N$, i.e.,
$$\mathbf{w}(n) = [w_0(n), w_1(n), \ldots, w_{N-1}(n)]^T$$
and the input vector $\mathbf{x}(n)$ should have the same length as $\mathbf{w}(n)$:
$$\mathbf{x}(n) = [x(n), x(n-1), \ldots, x(n-N+1)]^T$$
Therefore the output $y(n)=\mathbf{w}^T(n)\mathbf{x}(n)$ is a scalar, and the error signal $e(n)$ and disired signal $d(n)$ are also scalar. They represent the signals at discrete time $n$. You should update the filter weights sample by sample.


Which kernel produces sharper image kernel A = [ 0 -1 0 -1 4 -1   0 -1 0 ]
and B = [-1 0 -1  0 4 0 -1 0 -1 ]

i have found the frequency spectrums with freqz2 function of Matlab.

A's spectrum has the highest magnitudes in the corners (1,1),(-1,-1),(1,-1),(-1,1)
and B's in (1,0),(0,1),(-1,0),(0,-1).

I know that the center of the frequency spectrum has the low frequency information and the corners the high frequency information such as edges. It appears that B filter produces a sharper image although A has higher magnitude in the corners. Why is this happening ?

The sharpness of an image depends on acutance and resolution. The quantification of image sharpness seems image content dependent to me, since (i) sometimes the noise at the low SNR region may be enhanced with an high-pass kernel, and (ii) your A and B work on the edges with different orientations.

Take a very simple example:

T = [1 0 1;...
     0 1 0;...
     1 0 1];

conv(T,A,'same') gives you

   [4  -3  4;...
    -3  4  -3;...
    4  -3  4];

while conv(T,B,'same') results in:

   [3 0 3;...
    0 0 0;...
    3 0 3];

Looks like kernel A is better in sharpening the 45 degree edges. Yet with

  T = [0 1 0;...
       1 1 1;...
       0 1 0];

The performance of two kernels looks reversed (B appears to be better in horizontal and vertical edges enhancement).


Estimating scaling and translation between two consecutive temporal images I have two camera images taken from a camera mounted over a vehicle at time $t$ and $t+1$. I need to estimate the scaling and translation in  the images in spatial domain itself so that information from two images can be fused to get a better image. I have camera calibration parameters, vehicle dynamics parameters for the calculation.

Is it possible to estimate scaling efficiently? For translation I have used cross-correlation peak based approach, although it is not correct all the times.

I feel there is no the best method to solve this, but if you do not have some specific feature point to track, I suggest using optical flow a below:

compute optical flow of two frames,
find the dominate global optical vector, may be through computing average or median of optical vectors (you do not want small noisy motions)
Select at least 2 vectors parallel and similar to the global vector and fit a rotation*scale transformation matrix.

In case you have some feature points, instead of computing optical flow for all pixels simply use sparse optical flow OpenCV tutorial on Optical flow

This link might also be helpful.


Demodulating unknown signals - unmatched frequency sampling My goal is to make an approach to recognize modulation type of unknown signals.
I have energy detector which monitors spectrum and it provides information about center frequency and bandwidth of detected signal. After shifting desired signal to 0 frequency I make a resampling that new frequency sampling equals detected signal bandwidth. Problem is that I don't know exact sampling rate at which desired signal was generated at the transmitter. While analyzing the constellation plot after resampling there is no chance that convolution network will recognize modulation type.
Do you think that problem is unsolvable and this is waste of time?
Are there any techniques to overcome sampling frequency offset?

My goal is to make an approach to recognize modulation type of unknown signals.

Cool! Automatic modulation classification is still an exciting topic, after it was (one of) the foundations from which in the 1980s Software Defined Radio arose, e.g. [1].

Do you think that problem is unsolvable and this is waste of time?

Absolutely not! There's quite a few people and companies doing this, classically, commercially selling tools to automatically classify and demodulate signals [3], or SDKs with which you can build your own Machine Learning-based classificators [4].

Problem is that I don't know exact sampling rate at which desired signal was generated at the transmitter.

you don't really care about the sampling rate (can't tell that from the signal, ever; the transmitter contain a, say, 5/7-rate resampler that doesn't change the signal at all), but about the symbol rate.
What you hence are looking for is a symbol rate recovery algorithm, to know how often a symbol comes through, and a timing recovery algorithm, to tell you when exactly to look for the symbol. (The terminology I use here isn't that common, often "timing recovery" means both, getting the rate of symbols, and the exact temporal point.)
There's quite a few of these! The "classical" symbol rate recovery one for linear modulations (PSK, ASK, QAM, that kind of thing) is Oerder & Meyr [2], which exploits the fact that squaring the signal will introduce spectral components at multiples of the symbol rate. From there it's just a frequency estimation problem.
Problem is that it assumes a) sufficient SNR so that squaring doesn't make noise more powerful than your desired tones, and b) that the pulse shape after your receive filter adheres to the first Nyquist criterion. b) is a real problem, because you don't know the transmit filter, and systems usually only adhere to Nyquist I after the matched filter in the receiver (and that's matched to the unknown transmit filter).
There's some error considerations (and you can probably really work without exact knowledge of the transmit filter), but if you use a blind equalizer to bring your pre-channelized signal back to a white spectrum, you already have a starting point for your filter assumptions.
If you know anything about the transmitted signals (or dare to hypothesize and experiment), you can be much better. For example, if you know (or assume) the transmit signal is an OFDM variant, there's methods for that. For example, take this GNU Radio "Google Summer of Code" project [5], where a student designed a whole system to in general do something similar to what you describe you've done so far, and if it detects OFDM, to give you OFDM symbol length, bandwidth/rates, prefix dimension, and timing information.

[1] F. Jondral, "Automatic classification of high frequency signals," Signal Processing, vol. 9, no. 3, pp. 177–190, 1985; excerpt
[2] Oerder, M. ; Meyr, H.: Digital filter and square timing recovery. In: Communications, IEEE Transactions on 36 (1988), Mai, Nr. 5, S. 605 –612. http://dx.doi.org/10.1109/26.1476 . PDF – DOI 10.1109/26.1476. – ISSN 0090–6778
[3] PROCITEC: go2MONITOR, go2DECODE: Software for automated transmission classification and payload information extraction, https://www.procitec.de/en/products/
[4] deepsig: OmniSIG Sensor: detection and classification of RF emissions across large spectrum bandwidths in millisecond, https://www.deepsig.ai/omnisig
[5] Sebastian Müller: The Inspector. gr-inspector: A Signal Analysis Toolbox for GNU Radio, https://grinspector.wordpress.com/ . Demonstration video, Source code


High order IIR Comb Notch Filter with small Bandwidth - Low atenuation of harmonics Im filtering an EMG signal, more precisely the 50 Hz frequency and its harmonics. For that, I've constructed a comb notch filter (on filterdesigner in Matlab), with its fundamental frequency being 50 Hz, sampling frequency of 10kHz, and notches's bandwidth of 0.1 Hz. That means i have a filter order of 199.

Given the filter specifications, its magnitude response is bellow. I don't understand why the first two notches almost aren't attenuated. I've read that high quality factors (Associated with low bandwidth) are related with low attenuation rate, but i don't know in what extense it explains the low attenuation of the first two harmonics. Can you explain?

An order-200 IIR filter? You're playing with fire you know that? High-order IIR filter are prone to instability and/or misleading behavior even if they are stable.

If I understand correctly, you designed a 50-Hz notch filter, then a 100-Hz notch filter and so on?

Do you really need to filter out every harmonic? Do you know how high your harmonics are?

Solution #1 : Why not filter only the fundamental and the other major harmonics.  Usually even harmonics are not a problem and typically harmonics higher than 11 or 15 are small enough that you can simply ignore them. This would change your filter from an order 200 to a more manageable order 12 IIR filter.

Solution # 2 : Instead you could use a 200-tap moving-average that would allow you to extract the 50-Hz component and its harmonics. You could subtract the output of the moving-average from your signal and you would get rid of the 50-Hz component and its harmonics as long as they don't change too fast.


Relation between the matrix trace and the amplitude of each element Assume a diagonal matrix $\mathbf X$ whose size $N\times N$  and its diagonal elements are $0.5 + 0.5i$, and the vector $\mathbf p$ of size $N\times 1$ whose elements have similar amplitude.
I have noticed that
$$\operatorname{tr}\bigg\{\left(\mathbf{X^{-1}\cdot p\cdot p^H\cdot \left(X^H\right)^{-1}}\right)\bigg\}\tag{1}$$
is related into the amplitude of $\mathbf p$. where $\operatorname{tr}$ is the trace operator.

For example if $\mathbf p = 1$;
$$\implies\operatorname{tr}\bigg\{\left(\mathbf{X^{-1}\cdot p\cdot p^H\cdot \left(X^H\right)^{-1}}\right)\bigg\} = 2N$$

On the other hand, if $\mathbf p=0.5$;
$$\implies\operatorname{tr}\bigg\{\left(\mathbf{X^{-1}\cdot p\cdot p^H\cdot \left(X^H\right)^{-1}}\right)\bigg\} = \frac N2$$

And so on.

My question, is there a general formula expressing the relationship between the  equation $(1)$  and the amplitude of $\mathbf p$ ?

The trace of a matrix just scales with the scaling of the matrix:
$$\textrm{tr}\{c\mathbf{A}\}=c\,\textrm{tr}\{\mathbf{A}\}\tag{1}$$
So if you scale the vector $\mathbf{p}$, the trace scales with the squared magnitude of the scaling constant. With $\mathbf{p}=c\mathbf{p}_0$ you have
$$\begin{align}\textrm{tr}\{\mathbf{Ap}\mathbf{p}^H\mathbf{A}^H\}&=\textrm{tr}\{\mathbf{A}c\mathbf{p}_0c^*\mathbf{p}_0^H\mathbf{A}^H\}\\&=\textrm{tr}\{|c|^2\mathbf{A}\mathbf{p}_0\mathbf{p}_0^H\mathbf{A}^H\}\\&=|c|^2\textrm{tr}\{\mathbf{A}\mathbf{p}_0\mathbf{p}_0^H\mathbf{A}^H\}\end{align}$$


Implementations of Low-Pass Filter I should say which of the filters from the following image implements an ideal low-pass filter with cutoff frequency $\frac{\pi}{4}$. In the picture, Sa is a zero insertion filter, which inserts a zero sample after each input sample and Sb is a decimator, which extracts every second sample of its input. The first filter is fine and the second one does not represent the desired low-pass (that's the answer on my book).

However, I've not been able to explain why is the second filter wrong.
I know that up-sampling compresses the frequency content of a signal and that down-sampling expands it. So for the first filter, we first compress the input band by 2 and conclude that if the original limits were less than or equal to $\frac{\pi}{4}$, nothing gets filtered out (because $\frac{\pi}{4}$ becomes $\frac{\pi}{8}$ after Sa). And then, Sb recovers the original spectrum by multiplying by 2.

I analysed the second filter in a similar way, just noticing that Sb was switched with Sa thus making a multiplication by 2 first. The first Frequency Response would not filter anything from inputs with frequency lower than $\frac{\pi}{4}$ (which becomes $\frac{\pi}{2}$ after Sb). Then, after Sa, we have a division by 2 which restores the original spectrum and anything less than $\frac{\pi}{4}$ can proceed to the output.

What am I doing wrong?

If the input sequence is band limited to $\pi/2$, the second alternative is also OK (as you realized with your argument). However, if it not band limited you will get aliasing, so it is not a feasible solution in the general case.


MP3 Filterbank + MDCT: Why? What is the rationale for the two-step process the MP3 format performs, first decomposing the input into 32 subbands of 6/18 samples each then performing MDCT on each subband individually?
Why not simply perform a single MDCT?
Is this simply due to the non-power-of-two block size MP3 uses (no idea if there is a fast DCT algorithm for that size off the top of my head), or is there some more theoretical reason?

MP3 is a lossy compression which exploits -psychoacoustic masking.

Basically, the human ear has a perception limitations in identifying sounds due to active processing as well as due to non-linear frequency response.
Also note, that sound pressure level - the scale that corresponds to perceived loudness is actually logarithmic not linear.

OK - so this only indicates that we should "chop -off" frequencies that are not perceptible to human ear and keep those for which human ear has a high response. this simply can be achieved in 2 ways - either simply apply a filter to the sound and then encode it, or convert the entire signal into frequency domain and drop those frequencies that are on the edge.

NOW: MP3 does something much more than this.

Refer to the image here (Credit wikipedia):

What we see here, is that in the presence of frequency 250Hz, the 175Hz sound is not-audible at all to the human being unless that frequency also has dB level above a certain threshold mask. This frequency is called masking frequency. This is not a dominance of high frequency or low or otherwise. This is simply a function of which frequency produces sound at higher dB (at that time). the threshold of the audibility is also not static but depends on the dB level of the masking frequency. So in order to identify which all frequencies will go un-audible in the presence of other masking frequency one needs to first analyze the signal.

In actual effect the frequency masking effect is continuous in nature - but in practical codec - this is applied over bands.

Hence, we first divide the signal in frequency bands and then identify if any of the frequency bands will go un-audible in the presence of frequency bands which are louder. This is impossible to achieve through linear filtering of any kind! Hence, the filter-bank stage is employed.

UPDATE:

The different implementations of filter banks:
MP3 uses 2 stage filter bank - PQMF which produces 32 sub-bands followed by MDCT filter-bank which produces 575 bands. This is called hybrid filterbank. In general, MDCT can achieve higher resolution of filter bank which - in general PQMF can't achieve. Hence, Advanced Audio Coded (codec usually found in m4a and mp4 files) actually removes the redundant PQMF stage and uses MDCT event at a higher resolution. See Ref 3.

However, so the question is why is MP3 stuck to 2 stage filter bank?
As it turns out, the reason is more to do with compatibility with MPEG Layer I and Layer II. As indicated in the document here MP3 AND AAC EXPLAINED Ref 4
right from Fraunhofer, it says the real reason why MP3 remain stuck with 2 stage filter bank is due to backward compatibility with MPG 1 Layer I & II. From the paper (3.3.1):

  The polyphase filterbank has the purpose of making Layer-3 more
  similar to Layer-1 and Layer-2.

Ref [1]: Wikipedia - Psychoacoustics
Ref [2]: Theory behind MP3
Ref [3]: Modified Discrete Cosine Transform-Its Implications for Audio Coding and Error Concealment
Ref [4]: MP3 AND AAC EXPLAINED also available here.


Removing noise from PCM file I want to implement basic noise reduction. I'm recording audio using a microphone mounted on a drone. I would like to filter the audio file from drone noise.
What's the best way to do it ?

My first idea is to first record sound that is emitted by a drone only and convert it using fourier transform to see which frequency has the most energy.
Then I could do the same thing with audio file with drone noise.

My question is how to cancel this drone noise ?

This is the spectrum of the sound emitted by propellers:

Most of its power is in the low frequencies.

Edit.
Community noise, including aircraft and motor vehicles, typically ranges between 50 Hz and 5,000 Hz. So basically if I filter noise from propellers I might also filter sound from vehicles?

It really depends on the purpose of your project. If you want to do it from scratches, it will not be easy. However, if you want to do it for just one file, and the goal is to suppress the noise of the drone on that file, I recommend you to download free processing tool: audacity: http://www.audacityteam.org/download/, and use their feature where you can pick a period with a typical noise, analyze it and remove so analyzed noise from the whole file. I think that you may find their source code as well for free, which may include the implementation of the feature as well.
However, if you want to do it from scratches, there is a whole nice theory how to do it, but it is not easy.


Clarifying some notation in Bishop book In Bishop book, page 4, section 1.1, there's a notation I don't seem to understand what's meant by it. The whole paragraph, with which the section begins, is:

  We begin by introducing a simple regression problem, which we shall use as a running
  example throughout this chapter to motivate a number of key concepts. Suppose
  we observe a real-valued input variable x and we wish to use this observation to
  predict the value of a real-valued target variable t. For the present purposes, it is instructive
  to consider an artificial example using synthetically generated data because
  we then know the precise process that generated the data for comparison against any
  learned model. The data for this example is generated from the function sin(2πx)
  with random noise included in the target values, as described in detail in Appendix A.
  Now suppose that we are given a training set comprising N observations of x,
  written $x \equiv (x_1, . . . , x_N)^T$, together with corresponding observations of the values
  of t, denoted $t \equiv (t_1, . . . , t_N)^T$. Figure 1.2 shows a plot of a training set comprising
  N = 10 data points. The input data set x in Figure 1.2 was generated by choosing
  values of xn, for n = 1, . . . , N, spaced uniformly in range [0, 1], and the target
  data set t was obtained by first computing the corresponding values of the function
  1.1. Example: Polynomial Curve Fitting 5
  sin(2πx) and then adding a small level of random noise having a Gaussian distribution
  (the Gaussian distribution is discussed in Section 1.2.4) to each such point in
  order to obtain the corresponding value tn. By generating data in this way, we are
  capturing a property of many real data sets, namely that they possess an underlying
  regularity, which we wish to learn, but that individual observations are corrupted by
  random noise. This noise might arise from intrinsically stochastic (i.e. random) processes
  such as radioactive decay but more typically is due to there being sources of
  variability that are themselves unobserved.

I'm not clear what the notation  $x \equiv (x_1, . . . , x_N)^T$ might mean here, especially the so called variable T.
I guess it means it's the set of input vectors together with its corresponding set of target vectors (t). If I am right, so is this an example of  supervised learning? Still, what's T?

$x$ is a vector, composed of the individual elements $x_1$, $x_2$, etc.  The "T" is for transpose.  It is not a variable.


Duration of downsampled signal in analog frequency axis Suppose, I sampled a signal, $x(t)$ and let $x[n]$ be the sampled signal of $x(t)$. Let then I downsampled $x[n]$ by a factor of $M=3$. Let $x[n]$ after being downsampled by the factor of $M=3$ yields $y[n]$.
Now, my question is: if I plot the frequency spectrum of $y[n]$ along y-axis with frequency (in Hertz) along x-axis, then, will I get only $M=3$ copies of frequency scaled version of the original spectrum(spectrum of the original signal, $x(t)$)  or I will get infinite copies of the original spectrum ?

If you plot spectrum of y[n] along y-axis and frequency(in Hertz or in any other suitable  unit) along x-axis, there will be infinite copies of the frequency scaled(by a factor of M) version of the original spectrum(i.e, spectrum of original signal x(t)) each around kFs, where k is an integer and Fs = sampling frequency.


Deconvolution of shifted gaussian function in the frequency range I have a signal defined as
$$A(t)\cdot\exp\left(-i\omega_0t\right)$$
with $A$ the envelope function and $\omega_0$ the carrier frequency. I would like to transfer this signal into the fourier space and back, but need a large amount of points in time to accurately resolve the carrier frequency (the larger the time window, the more points are necessary). To circumvent that issue I intend to use a different strategy by using
$$F(A(t)\cdot\exp\left(-i\omega_0t\right))=\widehat{A}(\omega)\star\delta\left(\omega - \omega_0\right)$$
After the amount of points necessary to resolve $A$ accurately is significantly smaller than to resolve $\omega_0$ that significantly speeds up my calculation. Nevertheless, to apply the inverse Fourier transformation with the original amount of points used for the transform from $A(t)\rightarrow\widehat{A}(\omega)$ I have to "deconvolve" this function again by convolving it with
$$\delta(\omega + \omega_0)$$
but would require having negative frequencies. Is that even possible? If not, are there other ways to "shift" the function back, such that I can use fewer points for my FFT?

It is naive, but you didn't tell why you don't use $s(t) = A(t) \cdot \exp(-i \omega_0 t)$, and calculate $A(t) = s(t) \cdot \exp(i \omega_9 t)$...
I named your signal $s(t)$ here.
Other not so naive things.
You mention number of points, so I assume you are discretizing the signal. And your Fourier transform is a finite length DFT. In that case the Fourier transform will have alias, if you sample rate is enough to capture all the information in $A(t)$ then you can simply calculate $S(\omega) = DFT(s(k T_s))$, and the $A(\omega) = S(\omega + \omega_0~ \textrm{mod}~ F_s)$. Where $T_s$ is the sampling period, and $F_s$ is the sampling frequency, thus $s(k T_s)$ is the k-th sample of your signal.


IFFT of Microphone Array Directivity Response I am implementing a microphone array simulation in my thesis.
My prof and many documents indicates that microphone array behaves like filters in DSP but instead time-frequency domain there exists sensivity-directivity domains.
My microphone array directivity responce is a sinc like shape which can be seen below(x axis is angle from -90 to 90 ):

When I took the ifft of such directivity responce I am expecting to get a rectangle.
But instead ifft result is :

My microphone array composes of 41 apartures and distance between apartures is 10 cm.

My first question after taking ifft I am having 180 bins. How should I assign 180 values to my 41 apartures. Should I just take 41 bins in the middle?

My second question since microphone array directivity responce is a sinc like shape. I were expecting a rectangle like result after ifft. Why I am not getting rectangle?

Your plot shows the magnitude of the array gain vs angle of arrival for the broadside beampattern. Before you can take meaning IFFT of this sequence, you need to follow the following steps.

The data in the plot is the magnitude of the array gain, which is stripped of all phase information. For IFFT to give desired result, we need to restore the phase information too. Fortunately, this can be easily accomplished by just changing the sign of the alternate lobes, assuming the underlying curve is of $\sin(Mx)/\sin(x)$ form.
The gain of a linear array looking at angle $\theta_0$ is proportional to $\dfrac{\sin\left( Ma(\cos(\theta)-\cos(\theta_0))\right) }{\sin\left( a(\cos(\theta)-\cos(\theta_0))\right) }$. So it makes more sense to take the IFFT of the array gain sequence when it is plotted against $\cos(\theta)$ as opposed to just $\theta$, or angle (as it is shown in your plot).

Also, assuming that the measurement has been made at the design frequency (i.e. sensor spacing $d = \lambda/2$), the result of the IFFT would provide the sequence of sensor weights. The number of FFT bins would only affect the zero-padding around the actual sensors.

I took the liberty of trying the above steps out on a set of data shown below. Note that it is very similar to your top plot.

The second subplot shown the result as we restore the sign by inverting the sign of the alternating lobes.

At this point, let theta and gain represent the x and y data in this plot (after sign correction).

Next, the data is resampled such that the x-axis is linearly spaced in $\cos(\theta)$ domain. Note that since $\theta$ ranges from 0 to 180, $\cos\theta$ would range from -1 to 1.

NFFT = 2^8;
xr = linspace(-1,1,NFFT);
yr = interp1(cosd(theta), gain, xr, 'linear', 'extrap');

The choice of NFFT is arbitrary, but it must be greater than the number of sensors.

The result of resampling is shown in the third subplot above. Notice how it looks more like a sinc function now.

Finally, we are ready for the IFFT.

wgt = fftshift(ifft(yr, NFFT)); % weight of sensors
% plot
figure;
plot((-NFFT/2:NFFT/2-1), abs(wgt), 'o-', 'linewidth', 2);
grid on;
xlim(20*[-1, 1]);
xlabel('Sensor location (\lambda/2)');
ylabel('Sensor weight');
title('IFFT result (after FFTSHIFT)');

Based on the above plot, it appears that there were 13 sensors in the experiment.


The DTFT of $\{1,1\}$ is $1+e^{-j\omega}$ but what is the DTFT of $\{1,-1\}$? So I know that the DTFT of $\{1,1\}$ is equivalent to $1+e^{j\omega }$.

But what is the DTFT of $\{1,-1\}$ equivalent to?  Is it equivalent to $1-e^{j\omega }$?

You should consider the DTFT pair:

  $$\delta[n] \xrightarrow{\text{DTFT}} 1$$

and the time-shifting property of the Fourier transform

  $$x[n-n_0] \xrightarrow{\text{DTFT}} X(e^{j\omega})e^{-j\omega n_0}$$

plus the linearity of the FT. That is,

  $$ax_1[n]+bx_2[n] \xrightarrow{\text{DTFT}} aX_1(e^{j\omega})+bX_2(e^{j\omega})$$

Now consider the discrete signal $\{1,1\}$ which can be represented by $x[n], n\in\mathbb{Z}$ as follows
$$x[0]=1\\x[1]=1\\x[n]=0, \forall n\in\mathbb{Z} \ \backslash \ \{0,1\}$$

and can be denoted by $x[n]=\delta[n]+\delta[n-1]=x_1[n]+x_2[n]$. The DTFT is $$1+(1)e^{-j\omega(1)}=1+e^{-j\omega}$$
Similarly, the signal $\{1,-1\}$ can be represented by
$$x[0]=1\\x[1]=-1\\x[n]=0, \forall n\in\mathbb{Z} \ \backslash \ \{0,1\}$$

which can be denoted by $x[n]=\delta[n]-\delta[n-1]=x_1[n]-x_2[n]$, and the DTFT is $$1-(1)e^{-j\omega(1)}=1-e^{-j\omega}$$


Doing audio processing with DFT on a microcontroller I have written an audio processing application in C++ which runs on a desktop windows computer.

My next step is to implement this on an embedded system using a microcontroller, preferably arm-based.

As my application uses fft libraries which does floating point calculations, do I need a microcontroller with a floating point unit, or is there still a way to do this on a microcontroller without floating point hardware?

What type of ARM CPU would be recommended for processing audio with FFT like this? Thanks in advance

Generally, Arm Cortex-A9 type processor can be used for Audio signal processing blocks for effective results. For your case, basic arm processor type can support common floating point operations without any external floating point unit.
But, in case of special floating point operations like saturation,rounding and truncation, you shall use some advanced arm processors like cortex-a9 which has specific floating point unit for efficient signal processing.


Calculate the Inverse DTFT of the DTFT Derivative in Terms of $ x \left[ n \right] $ Consider the signal $ x \left[ n \right] $ and its DTFT transform $ X \left( {e}^{j \omega} \right) $.
Assume $ X \left( {e}^{j \omega} \right) $ is differentiable.

What is the Inverse DTFT of:

$$ j \frac{d}{d \omega} X \left( {e}^{j \omega} \right) $$

Express the result in terms of $ x \left[ n \right] $.

This is pretty straight forward using the definition of the Discrete Time Fourier Transform (DTFT).

The definition of the DTFT:

$$ X \left( {e}^{j \omega} \right) = \sum_{m = -\infty}^{\infty} x \left[ m \right] {e}^{-j \omega m} $$

Differentiating with respect to $\omega$:

$$\begin{align*}
\frac{d}{d \omega} X \left( {e}^{j \omega} \right) & =  \sum_{m = -\infty}^{\infty} \frac{d}{d \omega} x \left[ m \right] {e}^{-j \omega m} \\
& = \sum_{m = -\infty}^{\infty} \left( -j m \right) x \left[ m \right] {e}^{-j \omega m} \\
& = \frac{1}{j} \sum_{m = -\infty}^{\infty} m x \left[ m \right] {e}^{-j \omega m} \\
\end{align*}$$

Now, just plug it into the Inverse DTFT:

$$\begin{align*}
\frac{1}{2 \pi} \int_{0}^{2 \pi} \frac{1}{j} \sum_{m = -\infty}^{\infty} m x \left[ m \right] {e}^{-j \omega m} {e}^{j \omega n} d \omega & = \frac{1}{j} \sum_{m = -\infty}^{\infty} m x \left[ m \right] \frac{1}{2 \pi} \int_{0}^{2 \pi} {e}^{-j \omega m} {e}^{j \omega n} d \omega \\
& = \frac{1}{j} \sum_{m = -\infty}^{\infty} m x \left[ m \right] \delta \left[ n - m \right] \\
& = \frac{n}{j} x \left[ n \right]
\end{align*}$$

This implies:

$$ j \frac{d}{d \omega} X \left( {e}^{j \omega} \right) = j \frac{n}{j} x \left[ n \right] = n x \left[ n \right] $$


Regarding BIQUAD IIR implementation in MATLAB I want to implement a high pass filter (order 4) in MATLAB using cascade form.

[b,a]=butter(4,0.0083,'high');% I am using butterworth filter
[sos,gain]=tf2sos(b,a);

Should I use [sos,gain] for implementing the same in cascade form or just [sos]. I read a post related to same and it says just [sos] should be used, however, I am not sure if that is the correct way.
P.S. I want to create a MATLAB function for IIR filter(don't want to use the standard filter function in MATLAB) which would get the coefficients b,a & input.

EDIT: the answer below assumes that you called tf2sos with [sos,gain] as output arguments, and that your question is about the necessity of making use of the output variable gain. Thanks to Hilmar for pointing out this potential misunderstanding. His answer is about the equivalence of two different calls to tf2sos, one with two output arguments, and the other with just one. Note that when called with just one output argument, the resulting sos is different from the one obtained from a call with two output arguments.

If you want the correct scaling, in the sense that high frequency signals are passed through the filter with a gain of $1$, then you need to use the variable gain. If you don't use it you'll get the same high pass filter, but with some scaling. For fixed-point implementations, appropriate scaling is very important to minimize overflow and other quantization artefacts. For a Matlab implementation I wouldn't worry too much about scaling and just implement the filter with the correct gain.


How to create a band limited signal of sine waves that is nearly uniform distributed I would like to create a signal that is nearly uniform distributed (as depicted in Figure 1). It must consist of some single sine waves (e. g. 400 Hz, 430 Hz, 500 Hz) or a random band limited signal (e g. 495 to 500 Hz as depicted in Figure 2).
I expiremented with some parameters but I didn't get a good result. Does anyone know an analytical mathematical description or a kind of "thumb rule" for creating such signals?
Additionally I would prefer a signal that is as harmonic (in time domain) as possible (one single frequency would be the best - but of course not possible ;-) ). The signal is faded in and out with an raised cosine, because it's repeated periodically by a signal generator and I want to avoid steps in my signal.

To answer the "why" question:
I have an nonlinear system which I want to correct by avaraging (see picture 3). Usually I would get the red point which differs from the linear "true" value represented by the blue dotted line. The magenta colored dots are createt by a sine-wave and the big blue dot is gained through the avaraging of those magenta dots. This works whith a sine-wave signal, but since a sine-wave is not uniform distributet it is not perfect and the big blue dot differs from the big green dot representing the true value.
I tried to combine three sinusoids with different amplitudes. The result was not to bad but its a trial and error way. I would prefer an analytical way if there is any.

This is a wierd thing you are trying to do. What i the purpose behind this signal generation? What particular application needs this?

I think of two ways, but maybe one is not correct.

1 - Generate Uniform White noise, and then filter it to your desired bandwith of interest. But that will surely affect the probability distribution of the amplitudes and make it more laplacian, like your graph shows.

2 - Generate Any signal you want with the desired narrow bandwith (Still not sure why you need this, but let's say you do). you could use method one for it. Then estimate the probabiliy density with any method, like periodogram like you did, or Kernel Density estimation with parzen windows, etc.
Then, design a non-linear transformation that you can apply sample to sample, in order to go from the probability density that you have, to a uniform. This non linear transofrmation is actually the "Inverse Probability distribution" that you have.

If you estimate your density in discrete intervals with periodogram, maybe you could look at the form of the graph and try to guess what functional form it has, then fit it with least squares to a continuous function, and then obtain the Probability Distribution from that, find the inverse, and apply that transformation to every sample. That should give you a signal which is uniformly distributed.

I am guessing now, but i think that your bandwith will not be preserved by this transformation. I have actually never studied the relationship between bandwith and probability density of the amplitudes. But try it and see what happens.

hope it helps.

Edit:
I think you are trying to do is Least Squares fitting. It is the most mathematically correct application of "averaging to get rid of interferences and unwanted noise", or of linearizing a non-linear systems. Take all your measurments, and fit them to a line. There are 5 types of Least Squares methods:

Least Squares with $L_{2}$-norm
Least Squares with $L_{1}$-norm
Least Squares with $L_{\infty}$-norm
Total Least Squares
Data Least Squares

Try some of them and see which one gets better result for your particular application.


How to get the scale bar of an FFT for a 2d image I've been analyzing some images from a transmission electron microscope, including their FFTs, and I'm not sure how to apply a scale bar to the FFT images. I have calibrations for the real space image, and I was thinking it might just be the inverse of that scale, but I really have a hard time visualizing reciprocal space and would like to be sure.

This seems like it would have an easy answer, but I haven't been able to find it and I would appreciate your help.

If you are having trouble visualizing the reciprocal space, perhaps an analogy to 1D temporal FFTs will help (personally I'm more familar with these, so an analogy helps me anyway)

In a 1D temporal FFT, the input is N samples, at a rate of Fs samples/second, for a total time of N/Fs seconds.  The output is N samples, distributed between 0 and Fs (or between -Fs/2 and +Fs/2 depending on your point of view) (units 1/second).  So each output sample has width Fs/N (again units 1/second).

In a 1D spatial FFT, all we do is replace seconds by meters and samples by pixels. The input is N pixels, and a rate of Fx pixels per meter, for a total width of N/Fx meters.  The output is N pixels, distributed between 0 and Fx (or between -Fx/2 to Fx/2) (units of 1 / meter).  So each pixel has width Fx/N (agains units 1/meter).

So an output scale bar of M pixels has a width of M*Fx/N in units of 1/meters.

To go from 1D spatial FFT to 2D spatial FFT, just need a sampling rate in each direction.  Presumably you have the same sample rate in each directin.


What is an Intuitive Explanation of the Phase of a Signal I understand that the meaning of the phase response of a system is simply how much the system delays a frequency component. However, I do not find an intuitive explanation for the phase of a signal. The representation the Fourier transform of a signal by polar coordinates, produces two components; the magnitude and the phase. The magnitude response is simply how much each frequency component is contributing to this signal. However, what is the meaning of the phase here? Can I know if a frequency component is "exist" in the signal by just looking into the phase of a signal?

I found an explanation here:

  determines how the sinusoids line up relative to one another to form a
  signal

But the did not help me.

Compare $y= sin(\omega t)$ to $y=sin(\omega t-\phi)$ as shown in the plot below. $\phi$ is the additional phase term in radians where $\omega$ represents the frequency in radians per second. Thus the phase term shifts a sinusoid along the horizontal axis. So at a given frequency this will result in a time delay for that frequency, although to be noted that a fixed time delay will cause a phase shift that is proportional to frequency as a linear variation (The phase will increase linearly as the frequency increases). Thus our interest in "linear-phase" filters as this will have a constant "group delay" meaning the group of all frequencies in a waveform will be delayed the same amount in time, and thus not cause destructive interference between the different frequencies in the group. (Review Fourier Series and decomposition of waveforms into individual frequencies and the resulting reconstruction for more insight on that).

More pertinent (and useful) is complex signal representation with $Ae^{j \phi}$ representing the magnitude A and phase $\phi$ of a signal. This can be expanded using Euler's identity which may help give further insight: $$Ae^{j\phi}=Acos(\phi)+jAsin(\phi) = I + jQ$$

From which we see that in order to completely describe a complex signal we require two real signals, this could be magnitude and phase, A, $\phi$ as in $Ae^{j\phi}$ or the real and imagninary components I, Q as in I+jQ.

Refer to this post which may help give further insight into the importance and utility of complex representation for signals: Frequency shifting of a quadrature mixed signal


How to find the color channel shifting of an RGB image that was supposed to be greyscale? I am using a recording device that seems to shift the colors horizontally and I would like to find the amount of unshifting I have to do on 2/3 of the channels in order to obtain an image with minimal color distortion.

You can see here that in the original image at least the blue channel was shifted back with ~1.0 pixel.

The question is how can I detect the optimal unshift values, I'm sure they are not integers.

It seems like conceptually you're nearly there:

Take one of the three colour channels as reference (for example: red)
Match the green against the red
Match the blue against the red

The problem then is the matching with sub-pixel accuracy. To get
meaningfull numbers is a bit tricky, since the interpolation errors
any method will make have an impact on the accuracy.

There are two useable image registration methods I can think of:

Lucas-Kanade image registration. Using linear interpolation might
not give you accurate enough results, so consider bicubic or other methods.
Neil Dodgson has one made a nice overview. It's important that different
sub-pixel shifts of the interpolation kernel have a similar frequency transfer.
For the cubic family, the approximating b-spline is far better than cattmull-rom, in this case.
First upscale the image, then do pixel-precise image registration using your favorite
technique. (cross-correlation in the Fourier domain should do).
This only works if upscaling is done carefully. Bilinear or bicubic will most likely not give you enough accuracy. I can think of three ways:

a. Yen-interpolation. See equation 11 of his paper. Really slow, but optimal. Although you seem to have some alias in the signal so the 'bandlimited' assumption Yen makes might not hold.

b. Calculate FFT of the image, zero pad the high frequencies, inverse FFT.

c. Non-linear upscaling. Since the edges are quite sharp, the image is not properly bandlimited, that might be the main limitation of the previous methods. Edge Dependent Directional Interpolation might be better in this case.

Once you have the sub-pixel translation between the two, the correction to the green and blue channel is already solved, no matter which method you choose.


Under what circumstances is a non-uniform quantizer preferable over a uniform quantizer? I am studying btech ECE 3rd year. I want answer that question. Please tell me.

assuming the number of bits in the quantizer output word is the same for both uniform and non-uniform quantizer, then the number of quantization levels is the same for both.  then it is just a matter of placement of the quantization levels.

the "noise" or error of quantization is the difference between the quantized signal (the output) and the unquantized input.  the error will be proportional to the spacing of the quantization levels.  levels that are more closely spaces will have smaller quantization error than where levels that are spaced further apart.

let's assume the number of bits and the number of quantization levels are constant.

a uniform quantizer has uniformly-spaced quantization levels throughout the range.  large signals have the same spacing of quantization levels as do smaller signals.  that means that the smaller signals will have a lower signal-to-noise ratio (SNR).

a non-uniform quantizer can be designed so that the quantization levels are spaced more closely for smaller amplitudes and spaced more far apart for larger amplitudes.  as a result the signal-to-noise ratio can be made constant for both small signals and for large signals.  the SNR of the non-uniform quantizer would be worse than the uniform for large signals, but would be better than the uniform quantizer for small signals.


How are these pole zero plots created Back around 5 BC, where C stands for the programming language, engineers would photograph their oscilloscope's display to record their test results.

This appears to be such a photograph, taken from this Wikipedia article on Elliptic filters.

The photo is described as follows:

Log of the absolute value of the gain of an 8th order elliptic filter in complex frequency space. The white spots are poles and the black spots are zeroes. There are a total of 16 poles and 8 double zeroes. What appears to be a single pole and zero near the transition region is actually four poles and two double zeroes as shown in the expanded view below. In this image, black corresponds to a gain of 0.0001 or less and white corresponds to a gain of 10 or more.

I have seen plots like this many times and always wondered how they were created. In other words, what methods and equipment were used to do this? They even show poles in the RHP, which doesn't make much sense.

You can follow the links to the Mathematica source code, copy-pasted here:

xp2[xi_] :=
  Module[{g, num, den}, g = Sqrt[4*xi^2 + (4*xi^2*(xi^2 - 1))^(2/3)];
   num = 2*xi^2*Sqrt[g];
   den = Sqrt[8*xi^2*(xi^2 + 1) + 12*g*xi^2 - g^3] - Sqrt[g^3];
   num/den];
xz2[xi_] := xi^2/xp2[xi];

t[xi_] := Sqrt[1 - 1/xi^2];

(*Use particular values for low-order functions*)
r1[xi_, x_] := x;
r2[xi_, x_] := ((t[xi] + 1)*x^2 - 1)/((t[xi] - 1)*x^2 + 1);
r3[xi_, x_] :=
  x*((1 - xp2[xi])*(x^2 - xz2[xi]))/((1 - xz2[xi])*(x^2 - xp2[xi]));

(*Use nesting property for higher-degree functions*)
r4[xi_, x_] := r2[r2[xi, xi], r2[xi, x]];
r8[xi_, x_] := r4[r2[xi, xi], r2[xi, x]];

ellgain[xi_, w_, w0_, ep_] := 1/Sqrt[1 + ep^2*r8[xi, w/w0]^2];

DensityPlot[
 w0 = 1;
 ep = 0.5;
 xi = 1.05;
 min = 0.0001;
 max = 10;

 Log[Abs[
   ellgain[xi, sig + I*w, w0*I, ep]
   ]],
 {sig, -4, 4},
 {w, -4, 4},
 PlotRange -> {Log[min], Log[max]},
 PlotPoints -> 100,
 ColorFunction -> GrayLevel,
 ClippingStyle -> {Black, White}
 ]


Why Fourier transform and Stockwell-transform retain the absolute phase information of one signal? Hello friends am studying the topic of signal processing and the Fourier transform and the s-transform and in most books as for example "Time-Frequency Signal Analysis and Processing. 2nd" of Boashash in section 5.11.2.3 says that these transforms retain the absolute phase information of the signal and I do not understand that it mean with this concept about retain the absolute phase information of the signal and why this is important?.

It simply means that for a signal:
$$
x(t) = a(t) e^{j\phi(t)}
$$
that the phase information $\phi(t)$ is completely preserved in the transform
domain.

Time-frequency representations that use a quadratic form do not necessarily have this property:
$$
x(t+\tau)x^*(t-\tau) = a(t+\tau)a(t-\tau)e^{j(\phi(t+\tau) - \phi(t-\tau))}
$$
so if $\phi$ is of the form:
$$
\phi(t) = f(t) + C
$$
where $C$ is a constant, then $C$ is not preserved in $x(t+\tau)x^*(t-\tau)$.

This is important because sometimes you are interested in the absolute phase of a signal, not just its relative phase.


How to derive stability of analogue and digital systems I'm trying to find out why digital systems are stable iff poles oh their transfer function are inside the unit circle (or are there any others conditions?).

I understand analogous systems are stable iff all poles of their transfer function are negative 'cause of inverse Laplace transform, where $\mathcal{L}^{-1} \left \{ \frac {1}{s+a} \right \} = e^{-at} \cdot \mathcal{1} (t) $. (by the $\mathcal{1} (t)$ i mean heaviside function). Is that right?

But how to derive conditions for stability of discrete system?

What others ways for derivating stability of discrete and continuous systems exist?

In the discrete-time case you can use the $\mathcal{Z}$-transform:

$$\mathcal{Z}\{x(n)\}=X(z)=\sum_{n=0}^{\infty}x(n)z^{-n}\tag{1}$$

Equation (1) is the unilateral $\mathcal{Z}$-transform which is useful for causal signals which are zero for $n<0$. Now consider the signal $x(n)=a^n$, $n\ge 0$:

$$X(z)=\sum_{n=0}^{\infty}a^nz^{-n}$$

This series converges for $|a|<1$ and in this case $X(z)$ is given by:

$$X(z)=\frac{1}{1-az^{-1}}=\frac{z}{z-a}\tag{2}$$

From (2) you can see that $a$ is the pole of $X(z)$, and for $x(n)$ to be stable (if seen as a system's impulse response), and for $X(z)$ to exist, $|a|<1$ must hold, i.e. the pole must be inside the unit circle of the complex $z$-plane. This example is analogous to the continuous-time example in your question.

Also take a look at the concept of BIBO-stability (bounded-input-bounded-output stability), which leads to the same conclusion: for continuous time signals/systems the poles must be in the strict left half of the $s$-plane, and for discrete-time signals/systems the poles must be strictly inside the unit circle of the $z$-plane. The condition on the time-domain signals is that they be absolutely integrable in the continuous case, and absolutely summable in the discrete-time case:

$$\int_{-\infty}^{\infty}|x(t)|dt <\infty\\
\sum_{n=-\infty}^{\infty}|x(n)| <\infty$$


Understanding Cb and Cr Components of YCbCr Color Space I am familiar with additive (RGB), substractive (CMYK), and HSV-like colorspaces, but an article I'm currently trying to understand operates on YCbCr color space for image segmentation / object definition.

I've spend most of my morning looking for something that would explain the YCbCr naturally, but I just don't get it. I got a nice, intuitive explanation of the general idea behind this color space here, and explanation of how it's used for image coding/compression from these guys (all on photo.SE).

The formulas for calculating YCbCr from RGB is readily accessible on wikipedia.

I got the motivation for this representation, I got that Y component contains the most important (to the human eye) grey-scale information about the image.

I got that Cb and Cr carry information about the colors, and that (because of human eye (in)sensibility), they can be compressed without a visible lost in quality. But, what does each of the chrominance components actually represent?

As the article authors mention that "chrominance information is paramount in the definition of objects" in their approach, and I can not fully understand what I'm reading with my current "Y is intensity, Cb and Cr carry color information somehow" level of understanding YCbCr.

I'm seeking for an answer along the lines "Cb is..., while Cr is..." or "if you imagine looking through/with XY, you're actually looking at Cb component...", or some other way that would help me understand information carried by each of the components separately,  not just that they, together, carry color information.

EDIT

Let me give examples of intuitive explanations for other color spaces of the type I'm looking for:

RGB: Like shining a colored flashlight on a black wall: If you shine with a blue flashlight, you see a blue reflection. If you add a red flashlight, it will show a magenta reflection, which is a mixture of blue and red.

CMYK: Like mixing watercolors, you "add to the colors the surface reflects", (i.e. subtracts color from the background) so if you mix a yellow one with a cyan one, if will reflect green and thus you will get a green color.

HSV: Little kids are attracted to highly saturated objects, not bright (value). The Hue component is what "gives the color", while low saturation means the color is "diluted" by white. Change in value makes the whole thing brighter or darker.

With this definitions, I've been able to get an intuitive feeling about what a color representation in each color space means, without memorizing charts for each of them.

YUV (or YCbCr) is like HSV, but in different coordinates. (The difference between YUV and YCbCr is marginal - mostly related to exact formulas).

The $V$ component is the same. $(S,H)$ can be thought of as polar coordinates, and $(U,V)$ as cartesian. $H$ is the angle and $S$ is the radius. A rough conversion would be:

$ U = S \cdot \cos(H) $

$ V = S \cdot \sin(H) $

You can see this link for more information.

Another thing to add to your intuition list:

Saturation is how pure the color is from spectral point of view. For example, a laser has a very narrow spectrum, which implies high saturation.


Filtering out negative frequencies Is it actually possible to filter out the negative frequencies for the baseband signal is we take into consideration that the filter's impulse response is symmetric around the zero frequency?

I'm not sure if you ask about symmetry of the impulse response or symmetry of the frequency response. All filters with a real impulse response have conjugate symmetric frequency responses and symmetrical magnitude frequency responses. A symmetrical magnitude frequency response would affect the magnitudes of both frequencies $\omega$ and $-\omega$ identically so it could not possibly remove or attenuate a negative frequency while keeping the positive frequency intact. So the impulse response has to be complex.

The discrete impulse response $b_k$ of an ideal do-nothing filter is:

$$b_k = \frac{\int_{-\pi}^{\pi}e^{i\omega k}d\omega}{2\pi} = \left\{\begin{array}{ll}1 &\text{ if }k = 0\\\frac{\sin(\pi k)}{\pi k}&\text{ if } k \ne 0\end{array}\right.\\ = [\dots, 0, 0, 0, 1, 0, 0, 0, \dots]$$

The middle value $1$ is at $b_0$. Changing where the integral starts gives the ideal negative frequency removal filter:

$$b_k = \frac{\int_{0}^{\pi}e^{i\omega k}d\omega}{2\pi} = \left\{\begin{array}{ll}1/2&\text{ if } k = 0\\\frac{\sin(\pi k)}{2 \pi k} + i\left(\frac{1 - \cos(\pi k)}{2\pi k}\right)&\text{ if }k \ne 0\end{array}\right.= [\dots, 0, -\frac{i}{5\pi}, 0, -\frac{i}{3\pi}, 0, -\frac{i}{\pi}, \frac{1}{2}, \frac{i}{\pi}, 0, \frac{i}{3\pi}, 0, \frac{i}{5\pi}, 0, \dots]$$

The middle value $1/2$ is at $b_0$. The real part of the impulse response is zero at other values of $k$. The imaginary part has odd symmetry with respect to $\operatorname{Im}(b_0) = 0$. It is plain to see from the integral that the frequency response is zero for frequency $\omega < 0$ and constant ($1$) for $\omega > 0$.


Synthesizing harmonic tones with Karplus Strong I'm building my own MIDI synth.

I have to create a decent sounding tone for every note on the piano keyboard (A-1 through C7, or MIDI 21 through 21+88=109)

I've been using a basic ring resonator, where I create a ring of the right size (for example, sampling at 44.1 kHz the ring for A440 a.k.a. concert pitch would be 44100/440 i.e. 100.25 ~ 100 samples), fill it with static (values between -1 and +1) and then walk through it, averaging consecutive values.

The problem is that A-1 rumbles for 30 seconds, whereas C7 dies almost instantly.

I never noticed this problem before as I was working with a much smaller range; < two octaves.

Attempting to balance out the decays by multiplying by a suitable decay constant, say X[n] = .995*(X[n-1] + X[n-2])/2, where I choose the constant carefully for each note, is also failing.

Problem here is that even with decay=1.0, a high note still disappears almost instantly.  And it makes sense thinking about it; if the ring is only 10 samples,  even in 1/100sec we have covered 440 samples i.e. 44 revolutions, which will have smeared everything pretty close to 0.

One method I tried is feeding energy into the ring; creating a buffer of static enveloped by an exponential decay, and feeding this buffer into the ring.   It isn't working very well.

What I found works better is linear interpolation, so instead of doing ( X[n-1] + X[n-2] ) / 2, I do maybe:

0.98 * X[n-1]  +  0.02 * X[n-2] for a high note,  and
0.6 * X[n-1]  +  0.4 * X[n-2]  for a note lower down

etc.

But now the question is: how to choose the interpolation factor depending on the note?

I need a much better model to work with.

I think I also need to implement fractional delay;

Samp[ N ] += k *  Samp[ N - L ],

and if L = 10.2 say, I would have to do:

Samp[ N ] += k * ( 0.8 * Samp[10]  +  0.2 * Samp[ 11 ] )

so that is a completely different linear interpolation.

So my main question is: How to mathematically determine the necessary interpolation factor for a given note, so that all notes have an equal decay?

Supplementary questions are:
 - do I need to use a decay constant in addition to adjusting the interpolation?
 - is my understanding of fractional delay sufficient to get a working result?

If anyone knows of any resource for simulating a harmonic tone that is reasonably accessible, please do put it in as an answer. I've looked through JOS's resources, but he writes at a very advanced level, and I have difficulty understanding the material.

π

PS I can't find an appropriate tag;  pitch, note, gen*, tone, synth*, none of these are available.

First, note that the behaviour in which high-pitched notes disappear faster than low-pitched notes is found in many musical instruments. You might better not totally compensate for that.

One way I got the behaviour I wanted out of Karplus-Strong was the following:

I add an additional exponential decay for the lowest notes - a "loss" factor so that x[n] = loss * (x[n] + x[n - 1]) / 2
For the highest notes, I randomly decide whether I will do a buffer update or not ; with an increasingly high probability for the "no update" event as pitch increases.

To avoid dealing with fractional delays - and to get any kind of pitch modulation I want (but also to simplify the code), I use a ring buffer of the nearest power of 2 (for example, 128 samples for A4) and then do a linearly interpolated read from it. One can see this as a wavetable oscillator whose lookup table is dynamically modified/refreshed as the phase increment counter moves through it.


Magnitude and phase of $-\delta[n]$? I was reading this document and it shows the computation of the magnitude and phase of $h[n]=-\delta[n]$.

We can get the DTFT as:
$$H(e^{j\omega})= -1$$

So the magnitude will be $1$, and according to my understanding, the phase is calculated as $$\tan^{-1}\left(\frac{\Im\{H\}}{\Re\{H\}}\right) = \tan^{-1}\left(\frac{-0}{1}\right)=0$$

But in the document, the phase is $\pi$. Where is my mistake here?

If the magnitude is $1$ and the phase were $0$, you would get

$$H(e^{j\omega})=1$$

which corresponds to $h[n]=\delta[n]$. So this is obviously wrong.

The fact that the phase of a negative real-valued number equals $\pm \pi$ can be easily seen by drawing the number as a vector in the complex plane and measuring its angle with the positive real axis. Analytically, the phase angle of a complex number must be computed using the atan2 function. The formula with the $\arctan$ function that you used only gives the correct result if the real part of the complex number is positive. Note that by dividing the imaginary part by the real part, sign information is lost: if the result is negative you don't know if the real part or the imaginary part was negative, and if the result is positive, real and imaginary part must have the same sign, but we don't know which. This is why the $\arctan$ function of the ratio of the imaginary part and the real part is generally not sufficient for determining the argument (apart from the fact that you can't compute the ratio if the real part is zero).

In sum, a negative real-valued number $a<0$ can be written as

$$a=|a|e^{\pm j\pi}$$

Also take a look at this related question.


Structure Tensor vs. Hessian Matrix Hello can someone explain the semantical meaning of the structure tensor and the Hessian matrix. I am aware of how it is calculated, but i find it difficult to comprehend what they describe in the image and how they can be used to detect corners (Harris) and detect Blobs (SURF).

Any help is greatly appreciated, thank you in advance.

The structure tensor consists of first derivatives of the image. If the first derivative is high in one particular direction (one large eigenvalue of the structure tensor), then you have an edge. It the first derivative is high in two directions (two large eigenvalues), then you have a corner.

The Hessian matrix consists of second derivatives. Think of the Laplacian. It is the trace of the Hessian, and it is a center-surround filter, which gives you a high absolute value for a bright blob surrounded by dark background or vice versa. That's why it is blob detector.


Showing processes as an Equation I have developed an algorithm that counts the number of times a particular block (within a 2D block) crosses zero. Here's an example:

B1 = {  1, 4,
        -1, -3,
        0.12, -2.3
     }

The answer would there for be 3. To do this, what I do is:

Find the signum value of each element within the block (will return "1", "-1", "0") respectively.
If the signum value returns -1 then count increments by 1.

This will repeat until the there is no blocks, however, will only produce 1 value per block.

I am looking for a way to put all this process into an equation so I can demonstrate this rather than having to explain the processes in written text everytime. Is this possible?

Your algorithm, as described, is not correct.  You are simply counting the number of negative numbers, not the number of zero crossings.  For the particular data set that you use, B1, the two numbers are the same, but in general they won't be.

The way to do it is to multiply each number by the number after it, and if that number is negative, then you have a zero crossing.

zeroCrossings = sum((data(1:end-1) .* data(2:end)) < 0);

You also need to add some logic to handle samples that are exactly zero, but hopefully you get the idea at this point.


Why do we put poles close to the zeros in notch filter? In the book I'm studying from "Applied digital signal processing - Manolakis" it says:

  The problem with FIR notch filters is that the bandwidth of the
  notches is large. A simple way to create sharper notches is to place
  the zeros on the unit circle and two complex conjugate poles at the
  same angle with the zeros, close to the zeros but inside the unit
  circle.

But why is that? How it changes the thing? I tested it in matlab and in fact I get a way better filter but I don't understand the reason.

Consider the frequency response for a simple FIR notch filter given as a single pole and single zero:
$$H(z) = \frac{(z-z_0)}{(z-z_p)}$$
When the poles and zeros are plotted on the z-plane, the unit circle represents the frequency axis. So if we sweep z over all possible values where the magnitude of z = 1 (the unit circle) the resulting response will be the frequency response of the filter, both in magnitude and phase. Thus since the unit circle is $e^{j\omega}$ for $\omega$ in the range of 0 to $2\pi$, the frequency response is expressed as:
$$H(e^{j\omega})$$
That said we can easily now answer your question graphically (and thereby intuitively) if we consider the numerator and denominator  of $H(e^{j\omega})$ as two vectors in the 2d space of the z-plane as the value for z makes its way around the unit circle.  This is shown in the graphic below for the simple DC notch filter with a single zero at $z=1$ (which is DC), showing the numerator and denominator vectors when $\omega$ is little less than $\pi/4$: ($z=e^{j\pi/4}$).
The pole at the origin is "trivial" as in all cases the magnitude will be 1, so has no effect on the overall magnitude $|H(e^{j\omega})|$, but it does modify the phase of this transfer function. Therefore the overall transfer function will start with a null at DC and then grow rather slowly as z moves its way around the unit circle (all due to the magnitude of the numerator only).
The plot shows the pole and zero on the z-plane to the left, and then the vectors for the numerator and denominator are repeated on the right. (When subtracting 2 points on a 2D plane to describe a vector, the second point becomes the origin for the vector described).

Now for comparison consider what happens when we move the pole from the origin to be closer to the zero for the same position of z as in the graphic below. Notice in this case how the denominator is much smaller and quickly approaches the magnitude of the numerator, thus causing the DC notch to be much tighter.

The overall responses for both cases are shown below. The closer we make the pole to the zero, the tighter the response, but also the more precision is needed in a fixed point design.

This was explained simplest by using a DC notch filter. I have extended this to arbitrary notch filters for any frequency using the exact technique described by the OP (complex conjugate poles close to the complex conjugate zeros with the same angle) at this post:
Transfer function of second order notch filter


Best beginners guide to sound/audio I'm trying to learn audio programming, but am struggling with the basics of audio.

Audio programming is new to me (I'm a C/C++ developer with over 30 years of experience), but I don't have any background in sound/audio or DSP.

I'm learning ad-hoc, through various websites. Things I need to know are sampling, envelope, sweeps, length counters, pulse waves, triangle waves, delta modulation, PCM, etc.

Is there some kind of definitive book that I should start with, given my background? I don't really need to understand advanced sound processing, just something thats beginner and somewhat easy to digest.

In general you can divide the audio books into the following categories:

Audio Engineering
Audio Programming
Audio Signal Processing and Coding
Digital Signal processing

For your exclusive requirement (that fall into cat-2) I highly recommend the following book: The Audio Programming Book_BOULANGER_MIT


Where does the delta function come from if we derive autocorrelation directly? I am reading a book "Creating Noise" written by Hollos & Hollos and have a question about the autocorrelation function of the Gaussain white noise when reading the following passage:

From the theorem that the autocorrelation and psd are Fourier transform pair and the fact that psd of Gaussian white noise is $\sigma^2$, it is obvious that the autocorrelation of Gaussian white noise has a delta function $\delta(\tau)$ as formulated in (6). So at $\tau=0$, what we have is a pulse of inifinite amplitude. But if we derive the autocorrelation directly, as does the book in the passage, we only have $R(0)=\sigma^2$ which is a finite number, instead of an infinite pulse. I believe the autocorrelation result having a delata function is correct, but the direct derivation of it looks correct too. Being not able to find out where the mistake is, I am curious where the delta function $\delta(\tau)$ comes from if we derive the autocorrelation $R(\tau)$ directly from its definition. I wish you can give me some help with it. Thanks a lot.

The book is indeed not consistent. Continuous-time (zero-mean) white Gaussian noise has an autocorrelation function

$$R(\tau)=\sigma^2\delta(\tau)\tag{1}$$

and a constant power spectral density (PSD)

$$S(\omega)=\mathcal{F}\{R(\tau)\}=\sigma^2\tag{2}$$

The power (or variance) of such a process is

$$R(0)=\int_{-\infty}^{\infty}S(\omega)d\omega=\infty\tag{3}$$

So in the book the variable $\sigma^2$ denotes two different things: the finite value of the (constant) PSD, and the (infinite) variance of the process. But it can't be both at the same time. So the equation $R(0)=\sigma^2$ doesn't make sense, if at the same time we have $R(\tau)=\sigma^2\delta(\tau)$.

EDIT: In reaction to your comment concerning the definition of the auto-correlation function: the definition given in the book is correct (for real-valued wide-sense stationary processes):

$$R(\tau)=E[x(t)x(t-\tau)]\tag{4}$$

Since we know that $x(t)$ is white and zero-mean, we know that

$$R(\tau)=E[x(t)x(t-\tau)]=0,\quad\tau\neq 0\tag{5}$$

For $\tau=0$ we get

$$R(0)=E[x^2(t)]\tag{6}$$

which is the power of $x(t)$, and, since $E[x(t)]=0$, it also equals its variance. And now comes the error in the book: the authors use $\sigma^2$ to denote the variance: $$R(0)=\sigma^2\tag{7}$$ But then they claim

  The autocorrelation function of white noise must therefore be a delta function. $$R(\tau)=\sigma^2\delta(\tau)\tag{8}$$

But this is inconsistent with Eq. $(7)$ (of this answer). Note that both $(7)$ and $(8)$ are correct, but together they are wrong because they use the same symbol $\sigma^2$ with a different meaning. In $(7)$, $\sigma^2$ is the variance of $x(t)$, whereas in $(8)$, $\sigma^2$ is the value of the (constant) PDS of $x(t)$.

In the given example, $x(t)$ is white, so its power (variance) is infinite. This means that the constant $\sigma^2$ in $(7)$ cannot be finite. On the other hand, the constant $\sigma^2$ in ($8$) is finite, because it is simply the value of the constant PSD of $x(t)$.

The fact that the auto-correlation function of a zero-mean white wide-sense stationary process must be a weighted delta impulse at $\tau=0$ is obvious from the requirement that the PSD must be constant. The auto-correlation is the inverse Fourier transform of the PSD, and the inverse Fourier transform of a constant is a Dirac delta impulse.


Optimization of Filter Coefficients to Given Word Length It seems to me that this is a difficult mathematical problem where research is still carried on:
Is there an efficient way to find fixed-point coefficients whose zeros are closest to those calculated from the given floating-point coefficients?

It is definitely a difficult mathematical problem (integer programming), but I don't see much research going on in the DSP community. I do not know about the maths community. There were some important papers in the early 1980's on the optimum design of FIR filters with quantized coefficients:

D. M. Kodek, "Design of Optimal Finite Wordlength FIR Digital Filters Using
Integer Programming Techniques," IEEE Trans. ASSP 28, 304-308 (1980).
D. M. Kodek and K. Steiglitz, "Comparison of Optimal and Local Search Methods
for Designing Finite Worldlength FIR Digital Filters," IEEE Trans. Circuits Systems 28, 28-32 (1981).

But as far as I know none of the published methods developed into a standard design method for filters with quantized coefficients, unlike for the design of filters with infinite precision coefficients (e.g., Parks McClellan).

I think there are two reasons for that: first, optimum methods use integer programming techniques, which are computationally inefficient, and which will only compute a local optimum that can be arbitrarily bad. These problems become aggravated for large filter lengths. Second, there is not so much need anymore for the optimum design of filters with quantized coefficients. This has to do with the fact that more and more systems are implemented in floating point, and - if fixed-point is needed - the word length is usually large enough such that suboptimal heuristic design methods are sufficient. Such methods iteratively design the filter assuming infinite precision, quantize the coefficients, check the filter properties and, if necessary, adapt the specs (including the filter length, error weighting, etc.), re-design the filter according to the new specifications, and so on, until the requirements are met.

For IIR filters, the optimum design for quantized filter coefficients is even more difficult than for FIR filters due to the rational transfer function. I haven't seen any useful optimization methods in the DSP literature, only heuristic methods such as the one described above. Note that different filter structures will result in different possible pole locations in the complex plane after coefficient quantization. The possible pole locations for the direct-form structure are relatively sparse at small pole angles, i.e., for low cut-off frequencies, whereas the pole locations of the coupled-form structure are evenly distributed on a rectangular grid.


Upsampling and comb function Could somebody explain how upsampling works with comb function. How comb function looks for upsampling process?
(if it possible I would appreciate if answer would be related to DFT)

If with "Comb Function" you are referring to the periodicity of the spectrum of a sampled signal, first see Intuition for sidelobes in FFT for further details on the unfolded view of a digital spectrum (specifically how it repeats over the range of 0 to $F_s$ where $F_s$ is the sampling rate.)

(If you are referring to the use of "Comb Filters" then see this explanation of the CIC Decimator: Fast Integer 8 Hz 2nd Order LP for Microcontroller, and I can elaborate upon request on the CIC Interpolator which has a similar structure).

Upsampling and interpolation in general can be viewed as diagrammed in the figure below (once understanding the unfolded view of the digital spectrum from the top link above).

The frequency is upsampled by a factor of M by inserting M-1 zeros. In the case of the example in the figure, 3 zeros are inserted between every digital sample to increase the sampling rate by a factor of 4.

What happens when we do the zero insert is we get our original spectrum intact, along with new frequency components (replicas of our original spectrum) that now exist in our primary digital frequency span of $0$ to $F_s$. These are the exact same components that existed adjacent to all multiples of the sampling rate as described above in the unfolded view of the digital spectrum. If we can create a digital filter (Interpolation FIR or IFIR) that will pass the signal of interest with no distortion, and completely eliminate the unwanted components that would otherwise be distortion, as shown in the figure below, we will have a perfect interpolation between the original samples (the filter will grow the zero values to the perfect interpolation value needed, by nature of having all distortion frequencies eliminated while not distorting our spectrum of interest). Of course such an ideal filter is not realizable, but we can control and decide on how good of an interpolation we want based on the complexity (number of taps) in our filter design.

Also one may ask why the zero insert, why not do a first order hold (meaning use the current value as the value to insert instead of inserting zeros)?

To demonstrate the considerations with this, consider that a zero order hold is done mathematically in the time domain by convolving our waveform with a digital "pulse", while a zero insert is done by convolving our waveform with a digital "impulse".

Convolution in the time domain is multiplication in the frequency domain, so to explain using analog terms which the digital system will approximate- for the case of a first order hold, we are multiplying our desired frequency response by a sinc function (FT of a pulse) leading to passband droop, but for the zero insert, we are multiplying our desired frequency response by a constant (FT of an impulse), leaving our original passband undistorted.   The plots below demonstrate this for our intepolate by 4 example with the exact digital solution of what would occur.

Convolving with the digital equivalent of a sinc function is not a terrible way to go, and in fact its nulls will be right in the center of the frequency bands we want to suppress (very convenient). It is just important to understand the implications of passband droop and how to compensate for that if you take that approach (see how to make CIC compensation filter for a very simple 3 tap "inverse Sinc" compensation approach I have used). The CIC Interpolator is in fact such a "first-order hold" interpolation method that is very popular because of its simplicity.  Otherwise is absolute minimum distortion of your original waveform over its full usable bandwidth is important then zero-insert with a good filter design based on defined passband and stop band requirements is the way to go. Filter design is optimized by only rejecting the bands where the spectrum is replicated. Matlab's firls and firpm filter commands (Least Squared and Parks-McCLellan) both support creating multi-band filters which are perfect for this purpose. Further, if you go down the path of zero insert and optimum filter design, the end result can be implemented efficiently as a polyphase interpolator as shown in the figure below.  The process is to design the filter as described above, and then map the same coefficients into the polyphase structure as shown in the figure. The figure is descriptive of the mapping process in showing only 8 taps, but of course with only 8 taps in the original filter you cannot expect very good performance. An actual interpolate by 4 filter may have 64 taps for example, leading to 4 polyphase filters each with 16 taps, mapped the same way from the original base filter to the polyphase structure. The more taps the better the base filter can minimize pass band ripple and maximize the stop band rejections.

The reason for this mapping is very clear by following what occurs with the zeros that were inserted in the original filter- they mask several coefficients that are not needed on any particular given cycle (see figure below):


Is it possible the difference between OFDM-CP and OFDM-ZP about 10dB I'm building basic matlab code for ofdm, but when I change the CP OFDM into ZP-OFDM with sme parameters, I get the performance is worse by up to 10 dB!! Is that right ? or there is some thing wrong in the code !!

clear all; clc; clear;

N_sub = 64;   %Num of subcarriers
N_sym = 20;    %Number of symbols
Mod = 4;         %Modulation order
bi = log2(Mod);  %Number of bits
N_bits = N_sub*N_sym*log2(Mod);
CP = N_sub/4;  %CP length

%set channel parameters
h = randn(1,16);
h = [h zeros(1,N_sub-length(h))];
h = repmat(h,N_sym,1)';
ntap = size(h,1) - 1;

SNR_dB = 0:2:20;
errors = zeros(length(SNR_dB),1);

N_fram = 10^3;        %Number of frame
for snr = 1:length(SNR_dB)
    for i = 1 : N_fram
        bits = randi([0 1], 1 , N_bits);
        data_temp = bi2de(reshape(bits,N_bits/bi,bi));
        % Set modulation of data
        data_mod = (1/sqrt(2))*qammod(data_temp , Mod);

        data_fram = reshape(data_mod,N_sub,[]);    %% Serial to parallel conversion

        %FFt operation
        data_freq_total = sqrt(N_sub)*ifft(data_fram);

        %adding CP
        data_cp = [data_freq_total(end-CP+1:end,:); data_freq_total];  %Adding CP
%         data_cp = [zeros(CP,N_sym); data_freq_total];                  %Adding ZP

        %channel convolution
        for j_c = 1 : N_sym
            y(:,j_c) = conv(h(:,j_c),data_cp(:,j_c));                          % Linear convolution
        end
        y = awgn(y,SNR_dB(snr),'measured');                          % add noise

        %serial to parallel conversion
        y_s2p = reshape(y,N_sub+CP+ntap,[]);

        %remove cp
        y_cp2 = y_s2p(CP+1:end-ntap,:);

        %FFT operation
        data_freq_total_r = (1/sqrt(N_sub))*fft(y_cp2,N_sub);

        %Channel recovery
        y_fft_total = [];
        for j_r = 1 : N_sym
            y_fft_c = data_freq_total_r(:,j_r);
            y_fft_l = y_fft_c ./ fft(h(:,1),N_sub);
        y_fft_total = [y_fft_total y_fft_l];
        end

        %set demodulation
        data_rec = reshape(y_fft_total, N_sub*N_sym,1);
        data_demod = qamdemod(data_rec, Mod);
        bits_r = reshape(de2bi(data_demod,bi),1,[]);

        %calculate the errors
        errors_ofdm = sum(xor(bits,bits_r));
        errors(snr) = errors_ofdm + errors(snr);
    end
    snr
end
err_ber_fram_ofdm = errors/N_bits./N_fram;
figure (1)
semilogy(SNR_dB,err_ber_fram_ofdm)
grid on
legend('OFDM');

The difference I perform when modifying the CP-OFDM into ZP-OFDM is as below:

%adding CP
            data_cp = [data_freq_total(end-CP+1:end,:); data_freq_total];  %Adding CP
    %         data_cp = [zeros(CP,N_sym); data_freq_total];                  %Adding ZP

What's the issue of that degradation of performance ?

thank you

I tried using MMSE equalizer, it worked well. single tap equalizer is not always working with ZP as well as CP OFDM.


Inverse the Hanning Window I am using the following formula to calculate the Hann values for a set of points:
$$
w(n)=0.5\left(1-\cos\left(2\pi \frac nN\right)\right)
$$

To implement this in java, I am using the following code:

private static void applyHannWindow(final double[] points) {
    for (int i = 0; i < points.length; i++) {
        points[i] = 0.5 * (1.0 - Math.cos(2.0 * Math.PI * points[i] / points.length));
    }
}

To take the inverse, solve for $n$, I can re-use the $w(n)$ . Hence, the equation I come up with is this:
$$
n = \frac{1}{2\pi}\arccos\big(-2w(n) + 1\big)
$$

To implement this in Java, I am using the following code:

private static void reverseHannWindow(final double[] points) {
    for (int i = 0; i < points.length; i++) {
        points[i] = points.length * (Math.acos(-2.0 * points[i] + 1.0) / (2.0 * Math.PI));
    }
}

I wrote a simple test to make sure that after I apply the Hann window, and take it's inverse, that I get the same thing. However, I seem not to.

Given the original array:

final double[] data = new double[5];
data[0] = 200;
data[1] = 210;
data[2] = 215;
data[3] = 207;
data[4] = 203;

I get the following output:

### applying Hann window
0.0
0.0
0.0
0.9045084971874686
0.9045084971874787
### Reversing Hann window
0.0
0.0
0.0
1.9999999999999862
2.0000000000000138

Not sure what I am doing wrong here in applying the inverse to the Hann window?

I think you misinterpret the term inverse here.

Given a window sequence $w[n]$ of length $L$ with a region of support such that $ w[n] \not = 0 ,~~$ for $ 0 \leq n \leq L-1 $ , and a signal $x[n]$ you get the windowed signal $v[n] = w[n]x[n]$ , as each sample of the signal $x[n]$ is multipied by the corresponding sample of the window $w[n]$.

Inverting $v[n]$ means, given the windowed sequence $v[n]$ you obtain the original (unweighted) segment of $x[n]$ that resides in the region of suport range of the window $w[n]$ and that's given by: $$x[n] = \frac{v[n]}{w[n]} ~~,~~ w[n] \not = 0$$

The operation of recovering the original signal segment is as follows:
Let the original signal be the Matlab vector:

x = [1 2 3 0 5];   % original signal segment with 5 samples (L=5)

let the window be a Hamming type length of 5 samples:

w = hamming(5)';   % 5 sample Hamming window w = [0.08 0.54 1.0 0.54 0.08]

Hence the windowed signal is:

v = x.*w;          % v = [0.08 1.08 3.0 0.0 0.4] (windowed signal)

Now to recover x[n] from v[n] perform the following:

xrec = v./w;      % xrec = [1.0 2.0 3.0 0.0 5.0]; recovered

The operation $x_{rec}[n] = v[n]/w[n]$ is implemented per sample basis, such that $x_{rec}[0] = v[0]/w[0]$ or $x_{rec}[1] = v[1]/w[1]$ or in general
$$x_{rec}[n] = v[n]/w[n]$$

Now however the Hann window has the property that its initial and final samples are zero (in effect reducing the window length by 2) therefore an overlapping implementation is necessary to be able to recover $x[n]$ back from multiple $v[n]$ segments, as follows:

Let the original signal be the Matlab vector:

x = [1 2 3 4 5 4 3 2 1];   % original signal segment with 9 samples (L=9)

let's use a Hann window of length 5:

w = hann(5)';   % 5 sample Hann window w = [0.0 0.5 1.0 0.5 0.0]

First embed a zero to the first and last samples of x[n]:

xe = [0 x 0];   % xe = [0 1 2 3 4 5 4 3 2 1 0] (xe has 11 samples)

Now apply windowing into overlapping segments as such:

v1 = x[1:5].*w;          % v1 = [0.0 0.5 2.0 1.5 0.0] (first segment)
v2 = x[4:8].*w;          % v2 = [0.0 2.0 5.0 2.0 0.0] (second segment)
v3 = x[7:11].*w;         % v3 = [0.0 1.5 2.0 0.5 0.0] (third segment)

(notice the MATLAB indice operation x[i:j] to select portions of the input vector)

Now to recover $x[n]$ from $v_k[n]$ perform the following:
First note that each segment begins and ends in zero samples that you should discard, and you should consider only the remainig $L-2$ vetween those.

Then consider also the amount of overlap so that you don't recompute a previously computed input sample, then:

x1_rec = v1[2:4]./w[2:4];      % x1_rec = [1.0 2.0 3.0]; recovered first 3 samples of x[n]

x2_rec = v2[2:4]./w[2:4];      % x2_rec = [4.0 5.0 4.0]; recovered middle 3 samples of x[n]

x3_rec = v3[2:4]./w[2:4];      % x3_rec = [3.0 2.0 1.0]; recovered last 3 samples of x[n]


Why this DTFT plot from CTFT? I have a sample test with an answer but don't understand how they got to the answer:

  $x(t)$ has info only between $2 < |\omega| < 4$
  $X^F(\omega) = 0$
  for other frequencies.

  All noise was cleaned from the signal and its CTFT is defined:
  $X^F(\omega) = \begin{cases} |\omega|-2  & \text{for 2 < |$\omega$| < 4}\\ 0 & \text{otherwise} \end{cases}$

  The clean signal was sampled
  at $\frac 1 T = \frac {12}\pi$
  Draw the DTFT of $X^f(\theta) \text{ for} -3\pi < \theta < 3\pi$

  Use a clear coordinate system.

The answer has x-axis span of -9 to 9 and y-axis of 0 to 4.

How was this answer calculated?

I thought I would get a simple repeating copy of the original CTFT, but the original CTFT with $2 < |w| < 4$ is a line beginning at w=2 (on the x-axis) with the value zero, going up until the point  $w=4$ where it has the value of 2...

Please help!!  Thanks

When a continuous time signal is sampled in discrete intervals, the effect in the frequency domain is that the CTFT becomes periodic in frequency - that is, you end up with "copies" of the original CTFT, with each copy centered at frequencies equal to k*1/T (k = any integer).

So, when you sample at 12/pi you will get a copy of the original CTFT centered at w=0, w=12/pi, w=-12/pi, w=24/pi, w=-24/pi, and so on.

Since you have frequency content between 2<|w|<4 in the w=0 case, the w=k*12/pi cases will have frequency content between 2<|w-k*12/pi|<4.

Be careful - if frequency content from multiple copies have frequency content in the same places on the frequency axis, the effect is that they add together.  This is called "aliasing" and is usually an unwanted effect of sampling.

Here is a nice graph that shows the effect of sampling a continuous time signal on the FT.