Do FIR filter have any application in control theory? I know about IIR filter as discrete pendent to transfer functions in the Laplace domain. So it is actually quite simple to convert the function of the control system and finally receive the discrete output vector. (Bilinear Transform...)

But, how about FIR filter? Can we use them  to approximate functions of control systems.

FIR filters are fairly common in some areas of control theory. As they usually incur a lot of added phase/time-delay, they are not really usable in the feedback path of regular control systems, but they are useful when the added phase/time-delay is not affecting the system in an adverse way, or when the particular phase response and time-delay is desired.

Examples:

Feed-forward control. FIR filters are useful for producing filters that approximate arbitrary frequency responses, hence they can be used to shape a reference signal. A typical example is to use an FIR filter with the inverse frequency response of the plant -- trying to counteract the dynamics of the plant in order to get a desired output. Phase/time-delay is not interfering with the stability or performance since the computation can be done offline. FIR filters can often produce higher performance than IIR filters, especially where there are non-minimum phase zeros.
Learning systems (iterative learning control [ILC] and repetitive control [RC]). In these systems, again, FIR filters are often used to generate an inverse response filter of the plant. This is done to improve the performance of the learning, as this enables a higher gain to be used for a larger bandwidth. The time-delay is here desired (RC), or it does not matter, as the computations are done offline (ILC). Low-pass FIR filters are also used here for their linear-phase properties, as the linear-phase property makes the filter not interfere with the phase of other filters, and this improves performance.

FIR filters may be used to "clean" signals before they are processed by IIR. You would try to avoid noise from input the IIR, since it could lead to unstable states.


Periodicity of signals $$x(t) = \cos(2 \pi t) \cdot u(t)$$
$$y(t) = x(t) + x(-t)$$Is $y(t)$ periodic. If so, what is the $T$?
$$x(t) = \sin(2 \pi t) \cdot u(t)$$
$$y(t)= x(t) + x(-t)$$
Is $y(t)$ periodic?

Yes it is.
Y(t), X(t) have the same periodicity periode for the sin signal.
The cos signal is not periodic.

That is of course assuming that u(t) is the unit step function with u(0) =1.


DFT (FFT) of a Real Even Function Doesn't Yield Real Only DFT Signal There are a lot of similar questions already, but I can't find one answering this problem.

I used the following MATLAB code to perform the FFT on a Lorentzian input signal

x = linspace(-10,10,2^10);
y = 1./(x.^2 + 1);
F = fft(y);
IF = imag(F).^2;
power = trapz(IF); %integrating the imaginary part squared

Now, from everything I've read, the FFT of a real, even function is supposed to give $0$ imaginary part. But in this case, I get

power = 3.86

which doesn't seem to be just a small numerical inaccuracy. What's going on here?

Attached is the imaginary spectrum after fftshift.

When you work with DFT you need to remember one big assumption of the DFT - The samples are part of periodic signal.

For example, assume your signal in on the grid vX = [-3:3].
Let's say it is the simplest symmetric function: vY = vX .^ 2.
If you apply fft on it the result won't be real.

Why is that?
Because the DFT assumes the input is periodic.
Hence the input signal is something like [vY, vY, vY, ..., vY].

Moreover, for the DFT the indexing is [0, 1, ..., N - 1] which in this case is [0:6]. Now if we do periodic continuation, what will you get on the grid of [-7:6]?

figure();
plot([-7:6], [vY, vY]);

It won't be a be symmetric around 0 (Of course look at [-6, 6]).

Let's define vYY = vY(1:6).
If you apply fft(vYY) the result is pure real.
Why is that?
Because for the DFT it sits on the grid [0:5].
If you plot a periodic continuation of it:

figure();
plot([-6:5], [vYY, vYY]);

Look on [-5, 5], now it is symmetric.
Hence its DFT is pure real.

MATLAB Code

% The Grid
vX = [-3:3];
vY = vX .^ 2; %<! Even Function on the Grid

figure();
plot(vX, vY);
title('Figure 001');

vYDft = fft(vY);
norm(imag(vYDft), 'inf')

norm(imag(fft(vY(1:6))), 'inf')

figure();
plot([-7:6], [vY, vY])
title('Figure 002');

figure();
plot([-6:5], [vY(1:6), vY(1:6)]);
title('Figure 003');

Output:

  ans = 7.2678
  ans = 0

Figure 001

Figure 002

Figure 003

try

x = linspace(10*(-1024/1024), 10*(1023/1024), 2048);
y = 1./(x.^2 + 1);
y = fftshift(y);
Y = fft(y);
max(abs(imag(Y)))/mean(abs(Y))

i didn't run this, but i bet the number that appears at bottom is very very tiny.  like around 1e-15 .


FFT of sine wave not coming as expected i.e single point The cyan plot is a spectrum of 50 Hz, and the magenta one is a 50.1 Hz sine wave (having amplitude 0.7). Both are sampled at 1024 samples/s. I performed a 1024 point FFT to get this spectrum.

Why is only the 50Hz spectrum a single value? Why does the  50.1 Hz sine consist of other frequencies apart from 50.1 Hz; where do these new frequencies come from?

I did not do any non-linear processing on the 50.1 Hz signal! Also the 50.1 Hz appears to have smaller maximum amplitude, i.e. it's not 0.7, when in fact the sine wave I generated has an amplitude of 0.7.

Why is this?

Obtained by MATALB command fft();

Acually Matt's answer already gives one view on the problem here: the DFT is implicitly periodic in both time and frequency domain (see this question). From your parameters we can calculate that your observation period is 1 s. That means you observe 50 periods of a 50 Hz tone. Periodically extending that observation interval will always result in a seemless sine wave. If you take the 50.1 Hz tone, you're transforming 50.1 periods of an oscillation. Periodically extending that signal will result in phase jumps that cause additional spectral tributaries.

Another view is the frequency resolution that is defined as $f_\mathrm{s}/N_ \mathrm{DFT} = 1024 \text{Hz}/1024 = 1 \text{Hz}$. As 50 Hz is an integer multiple of 1 Hz, this tone can be reresented by a single DFT bin. As the result of the DFT is discrete values that are space 1 Hz apart, a 50.1 Hz cannot be represented by a single bin and the energy is smeared over the frequency domain. In your example this is the reason for the slightly asymmetric magenta spectrum and for the lower magnitude of the 50 Hz bin compared to the 0.7 magnitude of the 50 Hz spectrum.

Both effects described above contribute to the spectrum you're observing.

This occurs due to Spectral Leakage and Windowing. The ideal response i.e. impulse function is for continuous time sine wave. When you take DFT of a discrete sine wave in a digital computer, you are basically taking Fourier Transform of windowed and sampled sine and then sampling it in frequency domain. This causes the spectral leakage. Refer: http://w.astro.berkeley.edu/~jrg/ngst/fft/leakage.html


Instantaneous Frequency Estimation by Hilbert Transform - Theoretical Justification and Proof I would like to better understand why the instantaneous frequency estimation by Hilbert transformation works (and especially why it doesn't work / lead to precise results in many cases).

The motivation is to estimate signal $x(t)$ by decomposing it into an amplitude envelope $m(t)$ and phase of cosine $\omega_c (t)$ (or carrier waveform):

$$x(t) = m(t)  \cos\left(\omega_c (t)\right)$$

Now, assume that $x(t)$ indeed is a result from such a process.

Questions:

1) There are two "parameters" to be estimated for any $t$, as such some constraints are needed. What are the constraints regarding $m(t)$ and $\omega_c (t)$ that are selected when applying the Hilbert transform decomposition?

2) Is there a proof available somewhere that given the constraints, the estimation indeed finds the correct amplitude envelope and carrier (for continuous and also discrete case)?

See my earlier comments here: Meaning of Hilbert transform

Common fractal noise isn't an analytic signal (infinitely differentiable).  And a Hilbert transform re-creates the imaginary component of an analytic signal if you have the real component of the complex analytic signal (which one rarely has from real-world data).

However, sufficiently band-pass filtered data might be similar to a finite length segment of something that looks like an infinite length analytic signal (e.g. is from a source whose behavior can be modeled or estimated by a 2nd order linear differential equation).


basic hsb skin detection, neon illumination i hope this is the right place to ask. otherways sorry for my mistake and pleace advice me a better site.

i'm trying to implement a super simple skin detector using some range of hsb image. i'm using approach described here and here.

i'm try to use a video source from my webcam. if i use sun illumination it works quite well (not so good but quite good), but with neon light.. it is a mess. a lot of white region are detected and a lot of noise everywhere.

why?

i'm using the algorithm described in the second source:

convert image ho HSV color space
put white on the range 0 < H < 38
dilate filter
erode filter
blur filter

This actually might work better using a simple generative model in RGB rather than HSV.

Get a training image or several training images with some skin.
Manually select the skin pixels (e. g. by creating a binary mask)
Compute the mean and covariance of the skin tone in RGB (each should be 3-element vectors)
For an unknown pixel, compute its Mahalanobis distance from the mean, using the covariance.
Classify it as skin if the distance is less than a threshold.
Tune the threshold for best performance.

Edit:
I don't know if OpenCV has a function to compute covariance, but I can tell you how to do it yourself. Let's say you have $n$ RGB pixels.  You put them into an $n$ x 3 matrix, let's call it $P$.  Then compute $m$, which is the mean RGB vector by calculating the average of the columns.  $m$ will be a 1 x 3 vector. Subtract $m$ from every row of $P$ and call the resulting matrix $Q$.  Now to compute covariance, all you have to do is multiply $Q$ by the transpose of itself: $C = Q'Q$. Make sure that $C$ is 3 x 3.

Edit2:
The values you are getting seem to be too large.  To get the maximum covariance create the following matrix:

255 255 255
 0   0   0

and compute covariance of that.  You should get a matrix where every value is approximately
32513.  So make sure that your pixel values range from 0 to 255, and make sure you copy them into floats or doubles correctly.  The Mahalanobis distance is in the units of variance, so for that the numbers should be small.  Your threshold for skin classification should probably be less than 4.

import sys
import numpy
import cv2

cap = cv2.VideoCapture(0)
while(1):
    _, im = cap.read()

    im_ycrcb = cv2.cvtColor(im, cv2.COLOR_BGR2YCR_CB)

    skin_ycrcb_mint = numpy.array((0, 133, 77))
    skin_ycrcb_maxt = numpy.array((255, 173, 127))
    skin_ycrcb = cv2.inRange(im_ycrcb, skin_ycrcb_mint, skin_ycrcb_maxt)

    cv2.imshow("Second Image", skin_ycrcb) # Second image
    contours, _ = cv2.findContours(skin_ycrcb, cv2.RETR_EXTERNAL,
        cv2.CHAIN_APPROX_SIMPLE)
    for i, c in enumerate(contours):
        area = cv2.contourArea(c)
            if area > 1000:
                cv2.drawContours(im, contours, i, (255, 0, 0), 3)
    cv2.imshow("Final Image", im)         # Final image
    cv2.waitKey(1)


Modification of Accelerometer Data for Feature Extraction I have accelerometer [x,y,z] data recorded from various activities.  I am not looking at any filtering techniques (which I have seen within this forum) because I want to be able to measure the amplitude of the noise for each activity sample set using a global time-domain feature average loudness.

I am curious to know if I can center the signal around 0g, similarly to the X Data Minus Mean vector plot, but without the x-axis orientation i.e. the signal shifting up and down?

I have the x-axis data from one of my activities below:

I was hoping that by subtracting the mean from the x-axis data would yield this, but all it has done is shifted the entire signal down to around 0g.

The reason I ask is because I am unsure if computing a time-domain feature such as average loudness will work with the signal shifting up and down?

Average Loudness/Energy:

Your data is composed of a majoriy of low-activity regions, with higher-amplitiude  sparse/concentrated activities. The bad news is that most of the methods to achieve your goal are likely to somehow filter your data. The good news is that non-linear filtering could help you.

If you are not concerned with calculations and real-time yet, I would suggest first a centered median filter to extract a robust $0 g$ estimate, and remove it from the data. In other words: choose an appropriate left-right span of $K$. In each window around index $k$, computes the signal $\hat{s}(k) = s(k) - \textrm{median} [s_{k-K},\ldots,s_{k+K}]$.

A median is a non-linear filter minimizing the $ell_1$ norm (like the $ell_2$ norm minimization yields the standard mean). It is somewhat robust to outliers, and since yours are both above and below your thought $0 g$ estimate, a running median (computed around each sample) could extract this 0-level reference.

Looking at your plots, start with $K=50,100,200$ to see if this starts to provides you with what you are looking for. Then, more involded techniques could be considered.

The abrupt DC shifts appear like some trend. (Though off-topic, can I know what physical variation map to these DC shifts). As suggested previously, the problem can be circumvented, to certain extent, by using non-linear filtering. I can think about the following:

Carry out an Empirical Mode Decomposition (EMD, matlab code available here: http://perso.ens-lyon.fr/patrick.flandrin/emd.html). EMD decomposes the signal into few intrinsic mode functions (IMFs). When you add all the IMFs you get back the original signal. The key in the decomposition is the notion of data dependent filtering, and the last few IMFs capture the slow varying signal trend. Hence, I will suggest you do an EMD and discard the last few IMFs, and see the resulting summed up signal.
Differentiation: As suggested above but take care to choose your finite difference operation appropriately, else it may blow up the noise transients as well.
Though you are not willing to do this, I will still suggest, the average loudness here will refer to the high frequency spectrum energy content.This is implicitly what you will be obtaining from option 1 and 2.


A few overview questions about FFT. Relationships between input, buffer size, and expected output? I have a FFT set up and I've found that it responds very differently based on the buffer size and input. The output varies wildly, with magnitude of the peak frequency ranging from 500 to <1.

I have two tests in particular I'm running, the first being a sine wave at amplitude 1.0 with an LFO sweeping it's frequency at 0.1Hz. I've found that smaller buffers are more accurate in their ability to track the frequency. The magnitude of the frequency playing reaches 300 pretty easily.

The next test is with a mic, and it responds incredibly differently. If I whistle or hum, it will pick up that frequency. If I play music, rather than getting an increase in magnitudes of all frequencies, everything becomes weaker, including any whistling or humming. Changing the buffer size to be smaller helps, but at most iterations of the FFT on a buffer it results in very tiny magnitudes across the board.

In general, the things I'm finding are that smaller buffers are more likely to get a very strong magnitude of a frequency, but the more frequencies there are with a (magnitude > 1), the less likely that any of them are very strong (> 200).  But none of this really means anything if it's a product of a broken implementation, and I can't seem to find any specific numeric information on what my expectations should be.

So, given an audio input and an output of frequencies and magnitudes (calculated by $\sqrt{\frac{\text{real}^2 + \text{imag}^2}{\frac{\text{numSamples}}{2}}}$), my questions are as follows:

Given a sole, unmoving sine wave, with a relatively small buffer size what kind of magnitudes are to be expected on what frequencies?
Given a sole, unmoving sine wave, with a relatively large buffer size what kind of magnitudes are to be expected on what frequencies?
Given a sole, sweeping frequency sine wave, with a relatively small buffer size what kind of magnitudes are to be expected on what frequencies?
Given a sole, sweeping frequency sine wave, with a relatively large buffer size what kind of magnitudes are to be expected on what frequencies?
Given static, with a relatively small buffer, what kind of magnitudes are to be expected?
Given static, with a relatively large buffer, what kind of magnitudes are to be expected?

I know this is a lot of questions, but the results I'm getting from these tests are not what I was expecting, and I want to start by making sure it's not my expectations that were wrong.

So, I'm not an expert at all, but I feel that I've perhaps learned enough to be somewhat helpful to anyone who ends up here with little introduction to FFT.  If any of this is wrong, see the "P.S" at the bottom.  tl;dr I'd rather share than leave this question unanswered entirely.

A few things are important to note: Parseval's theorem, bin frequencies, and normalization:

Parseval's theorem is something to the ends that the energy of the signal will be the energy of FFT bins (bins being each frequency/magnitude pair).  So, what this means in practice (I think) is that a small buffer size will make the bin magnitude of a sine wave frequency vary based off the amount of high amplitude wave is in the buffer.  It won't be consistent.

Bin Frequencies are the frequency of each bin from your FFT.  They aren't going to be the exact numbers you are looking for without a bit of work, but they should be close enough to give you some info.

Normalization is something that helps with all this.  For my project, I found the max bin magnitude, averaged it with the max from the last cycle, and then divided everything by it.  The result is that my visualization is showing something much more like what my ears feel they would see.

Now for the Q's (assuming all signals are normalized to -1 < x < 1):

Magnitudes of ~200 +-100, in the bins near the signal wave freq
Same, but magnitudes of ~300 +-50.  The larger buffer limits cases of more or less high amplitude chunking (sin(x) at [0 to Pi/4] vs [Pi/4 to Pi/2])
Magnitudes of ~150 +-100, in the bins near the signal wave freq.  The sweeping spreads out the bins a little, but the small buffer size (relative to sweep speed) helps.
Magnitudes of ~50 +-100, in the bins touched by the signal wave freq.  If the sweep is larger or faster expect the magnitudes to be even more spread out.  Note: the more bins touched, the lower the magnitudes overall.
Magnitudes of ~5 +-5, in all bins which happen to be strongly represented in this small buffer chunk
Magnitudes of ~1 +-1 in all bins.

So, it seems there is probably some form of log function which would help make the perceived levels of all these sounds not so distant when graphed by magnitude, but, overall, using normalization has worked in my project to pick out the most noteworthy note tones.  (btw notes and tones are anagrams)

P.S.  If any of this is wrong, please do make an answer covering and correcting these topics.  I will not be marking this as the accepted answer for at least a month because I'm not experienced enough to be an authority on this subject.  That being said, this is what I think may be true based off the results I've been getting, and it's better than nothing for a newcommer.


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.

Separating the input signal into subbands allows to use different MDCT parameters for each band.

These parameters are determined by the psychoacoustic model, which determines, for each subband, the window length for the DCT and quantization steps for the resulting coefficients. See, for example, slide 9 of this summary.

Performing only one DCT for the whole signal would not allow the psychoacoustic model to adjust the admissible quantization noise and window in each subband, which is the basis for data rate compression.

Further compression is then performed by source coding (Huffman).


What is the exact meaning of unstable system in DSP? In physical systems I understand what is the meaning of stability or unstability. An operational amplifier for example, if working in positive feedback will either saturate or start osscilating (i.e will not have any stable state). thats clear to me.

But I am unable to understand what exactly we mean when we say an IIR filter (or any other digital system), for example, can become very unstable.

What exactly inside the Digital Signal Proccessor  happens, what happens to the output physically?
What exactly we mean by unstable system in this context?

Unstable typically means and unbounded output for a bounded input. In other words the output of, say , a filter can get infinitely large although the input is perfectly okay and of "normal" size.
A simple example would be the difference equation $y[n] =  x[n] + y[n-1]$. If we calculate the step response, i.e. $x[n] = u[n]$, we get y[0] = 1, y[1] = 2, y[2] = 3 ...
The output grows infinitely even though the input is a perfectly well behaved signal, bounded by 1.

When a system is unstable, the output of the system may be infinite even though the input to the system was finite. This cause a number of practical problems. For instance, a robot arm controller that is unstable may cause the robot to move dangerously. Also, systems that are unstable often incur a certain amount of physical damage, which can become costly. Nonetheless, many systems are inherently unstable - a fighter jet, for instance, or a rocket at liftoff, are examples of naturally unstable systems. Although we can design controllers that stabilize the system, it is first important to understand what stability is, how it is determined, and why it matters.


Vibration Sensor - Balancing Machine Im' working on a balancing machine. this machine has two vibration sensors. (pizzo) my question is about that why this machine use two Sensors? this machine has an rotary encoder to detect disk angle.

i trying to read these signal by ADC card and using FFT bandpass filtering to calculate unbalanced mass on disk. for this algorithm i just use one of Vibration sensors. should i use second sensor?

Of course, the machine could be provided with three, four or more sensors, if the goal of the control unit is perform a balancing of the rotation. Actually three sensors should be the proper selection if the rotation is unconstrained in a biaxial plane, and two sensors, if constrained in a uniaxial direction.


How to get a generic 3d Face model? I'm working on image processing that involves facial features.

Where can I get a generic 3D model of a human face? (any format is fine as long as I can manipulate the data using C or Matlab)

You might have to make your own from a 3d face database, which may not be easy to do. This website has a good listing of public databases. http://www.face-rec.org/databases/

A Open-Source software called MakeHuman might be what you are looking for.


How to verify the obtained SNR of the noisy signal after adding noise? I defined the SNR as follows:
$$\rm SNR = \frac{\text{Signal power density}}{\text{Noise power density}}$$
Given noise-free signal, I added the noise to this signal by using the formula above. Now, I am looking for a way to verify if the obtained SNR of the noisy signal is equal to the SNR value that I set when I generated the noise.
Thanks in advance!

Here is a common/straight-forward way using the discrete time domain samples. If you have the noise free signal, $x[n]$, and you created a noise signal, $w[n]$, then you can calculate the SNR by using the formula:
$$\hat{\text{SNR}}=\frac{\sum_n |x[n]|^2}{\sum_n |w[n]|^2}$$
Edit
An unstated assumption above is that $\mathbb{E}\big[x[n]\big]=\mathbb{E}\big[w[n] \big]=0$. In general, you'd use the variance not the sum of squared values, as to take care the possible non-zero means:
$$\hat{\text{SNR}}=\frac{\text{var}(x[n])}{\text{var}(w[n])}$$

On SNR wikipedia page, there is a section "Alternative definition".
SNR = Mean^2 / Standard_deviation^2


Bandwidth of a composite signal How to find bandwidth of a composite signal that can be decomposed into five sine waves with given frequencies and having all peak amplitudes same ?

Unless I'm missing something here it would just be the difference in frequency between the highest frequency sine wave and the lowest:

BW = Fmax - Fmin

The bandwidth in this case would be the maximum frequency of the five sine waves. There is no problem if you know the five frequencies a priori. However, if that information is not available, you can do a DFT to find out the frequency content and find the index upto which the frequency content is non zero. That would give you your bandwidth.


Sampling Theorem illustration Can someone please explain the illustration (figure 1.19) at the bottom of this image? It looks like there are four sampling points but I don't understand what the different curves represent. Actually except from the dashed curve I can't even distinguish the different curves from each other. Can someone please clarify this so that the illustration makes sense?

(Image source: Digital Signal Processing Principles Algorithms and Applications)

This plot depicts how to convert your digital signal back to the analog one, using $\mathrm{sinc}$ functions. The nice property of these functions used in this process, is that maximum of each function always occurs at minimums of the other, shifted function:

Now the process of the D/A conversion that is depicted on your plot is basically taking a $\mathrm{sinc}$ function, and multiplying it by one sample that is at the maximum of $\mathrm{sinc}$. You do it for each and every sample of your digital signal, and then sum all of the outputs. For example, for the above plot, the outcome would be following:

Now maybe some more complicated example of D/A conversion. On the plot below, you can see the actual samples (in blue) and "analog" $\mathrm{sinc}$ interpolated signal (in red).

Plotting each function separately you can observe what will become an analog signal later:

I hope that it is more intuitive for you now. In order to complete this example, here is some MATLAB code I used:

clc, clear all, close all
fs = 32;
dt = 1/fs;

% Generate a samples of the signal
ts = linspace(0, 1, 10);
s = 2*sin(2*pi*ts) + sin(2*pi*3*ts);

td = linspace(0, (length(s)-1)*dt,     length(s));
t  = linspace(0, (length(s)-1)*dt, 100*length(s));

sincM = repmat( t, length(td), 1 ) - repmat( td', 1, length(t) );
SS =  s* sinc( sincM/dt );

% Sampled signal
stem(td, s, 'b', 'linewidth', 2)
hold on
% Reconstructed signal
plot(t, SS, 'r', 'linewidth', 2)
grid on

legend({'sampled signal', 'interpolated signal'})

figure
stem(td, s, 'b', 'linewidth', 2)
hold on
grid on
% Each sinc function
basis = sinc(sincM/dt)';
sbasis = bsxfun(@times, basis, s);
plot(t, sbasis)

The explanation given in your text is, unfortunately, quite dreadful
and confuses several issues. We have a low-pass function $x(t)$ whose sample values $$\cdots, x((n-2)T), x((n-1)T), x(nT), x((n+1)T), \cdots$$
are given to us. Note the complete absence of any subscript $a$ on $x$.
We write $x_n$ for the $n$-th sample value $x(nT)$. We wish to reconstruct $x(t)$ from its samples, and our reconstruction
is $x_a(t)$ as (incorrectly) defined in (1.4.23): the correct
definition is
$$x_a(t) = \sum_{n=-\infty}^\infty x(nT)g(t-nT)
= \sum_{n=-\infty}^\infty x_n g(t-nT).\tag{1.4.23.corrected}$$
Note again the complete absence of any subscript $a$ in
the middle sum, and the right-hand sum emphasizes that the reconstruction
is a weighted sum of delayed versions of $g(t)$ with the weights
being the given sample values.

The function $x_a(t)$ as
defined in $(1.4.23.\text{corrected})$ has the property that if $t=mT$,
then
$$x_a(mT) = \sum_{n=-\infty}^\infty x_ng(mT-nT)
= \sum_{n=-\infty}^\infty x_n g((m-n)T) = x_m \tag{1}$$
because $g(t)$ is a Nyquist pulse and enjoys the property
that $g(\ell T) = 0$ if $\ell$ is a nonzero integer while
$g(0) = 1$. Thus, the reconstructed function $x_a(t)$ has
the property that at each sampling instant $mT$, the value of
$x_a(t)$, viz. $x_a(mT)$ equals $x_m$, the $m$-th sample
value given to us.  Thus, as far as samples are concerned
the two functions $x_a(t)$ and $x(t)$ have the same values
at the sampling instants $mT$. In other words,
$$x_a(t) = x(t) \quad \text{for} \quad t = mT, m \in \mathbb Z.\tag{2}$$
Note that $(2)$ holds regardless of what $T$ is and how it
relates to the maximum frequency in $x(t)$ etc.

What at time instants other than multiples of $T$? Is it
true that $x_a(t)$ as defined in $(1.4.23.\text{corrected})$ equals
$x(t)$ for all $t, -\infty < t < \infty$?
Well, that is where the fun statrts. If the Nyquist criterion is met, then $x_a(t)$ as given in $(1.4.23.\text{corrected})$ is precisely
the $x(t)$ that we wish to recover. I will not prove this statement
(which is called the Nyquist theorem or the Shannon-Whittaker
sinc interpolation formula etc) but that $x_a(t)$ as defined
in $(1.4.23.\text{corrected})$ is $x(t)$ is what the crux of
the matter is.  If the Nyquist criterion is
not met, then $(1.4.23.\text{corrected})$ gives a function $x_a(t)$
whose sample values
$x(t)$ at the sampling instants $mT$ (cf. the remarks surrounding
$(1)$ and $(2)$ above) but is not equal to $x(t)$
at other times. In fact, $x_a(t)$ is the result of aliasing
of $x(t)$ with respect to the sampling frequency etc, but that
is a matter for a different question.

With that as preliminaries, let us turn to the illustration that
puzzles you. The dashed curve is the result of the interpolation:
the function
$x_a(t)$ as defined in $(1.4.23.\text{corrected})$ and it is
the sum of infinitely many weighted and delayed sinc functions.
The other curves are (four of the)
infinitely many sinc functions $x_m g(t-mT)$ on the right side
of $(1.4.23.\text{corrected})$. At $t=nT$, the term $x_ng(t-nT)$
which occurs on the right side of $(1.4.23.\text{corrected})$
has a peak value of $$x_ng(nT-nT) = x_ng(0) = x_n$$ and it has
value $0$ at other sampling instants
$\cdots, (n-2)T, (n-1)T, (n+1)T, \cdots$ Similarly, the function
$x_{n-1}g(t-(n-1)T)$ has value $x_{n-1}$ at $t=(n-1)T$ and value
$0$ at all other sampling instants.  In short, if we evaluate the sum
$$\sum_{m=-\infty}^\infty x_m g(t-mT)$$
at $t=nT$, that sum of infinitely many terms simplifies
to $x(nT)$ because all
the other terms are $0$. No muss, no fuss, no worries about
convergence of an infinite series etc.  Remember that
all this holds without any
issues of Nyquist criteria etc. The hard part is the assertion
that $x_a(t)$ as defined in $(1.4.23.\text{corrected})$ equals
$x(t)$ for all $t$ if the Nyquist criterion is satisfied.

So why all this yelling and screaming about (1.4.23) being incorrect
and the explanation in the text being sheer nonsense?  Well,
(1.4.23) as shown in the OP's text has sample values of $x_a(t)$,
the reconstructed signal, on the right side. So, the definition
is circular, to say the least, and even if we ignore that, it
is an assertion of the Nyquist theorem or the Shannon-Whittaker
interpolation formula without any proof. The equality
(1.4.23) as given in the text holds at $t=mT, m \in \mathbb Z$;
the assertion that it holds for all $t$ is a proof by intimidation.


Detecting a filled glass object I've been sent here from this question in stackoverflow, please excuse me if the question comes too specific and it's not in the manners in here:)

The task is to find a glass with specific liquid in it. Let me show you the pictures and then describe what i'm trying to achieve and how i was trying to achieve so far in the description below the pictures.

The pictures: (seems i need at least 10 reputation to post pictures and links, so links will have to do :( otherwise you can look at the stack overflow question)

A detailed description:
I was trying to implement an algorithm that would detect a glass of a specific shape in opencv (glass may be transformed by a different camera shot angle/distance). There will be also other glasses of other shapes. The glass i'm searching for will also be filled with some colored liquid that will distinguish it from glasses containing other colors.

So far, i have tried using SIFT feature extractor to try to find some features in the glass and then match them with other photos with the glass in it.

This approach worked only in very specific conditions where i would have glass in a very specific position and the background would be similar to the learning images. The problem also is that the glass is a 3d object and i don't know how to extract features from that (maybe multiple photos from different angles linked-somehow?).

Now i don't know what other approach could i use. I have found some clues on this (here https://stackoverflow.com/questions/10168686/algorithm-improvement-for-coca-cola-can-shape-recognition#answer-10219338 ) but the links seem to be broken.

Another problem would be to detect different "levels of emptiness" in such glass, but i haven't even been able to find the glass itself properly.

What would be your recommendations on the approach in this task? Would it be better to use a different way to find the local 3d object feature? Or would it be better to use other approach altogether? I have heard about algorithms "learning" the object from a set of multiple photos but i have never seen this in practice.

Any advice would be really appreciated

The paper referenced in your link seems to be this one.

Of particular interest there is Table 1 (included below). The accuracy rates aren't great, though they are better than other approaches.

There has been much work done on this when it comes to facial recognition software. For example, if you notice on Facebook when tagging photos the location of the faces are boxed and suggested to you.

I have seen a lot of literature on facial recognition in images using neural networks and a quick Google search will undoubtedly turn up a vast amount of information on the subject. These networks take the pixels of the image as inputs. In your case, the way the opacity changes/light reflects off the glass may be good identifying features that the network will learn.

One problem may be the number of photos you have to use as training data and the pre-processing of these (i.e. identifying the faces yourself). If it is unfeasible to do this for enough images to train your network well enough then you will have to look for some shortcuts in the learning stage. This paper is relevant to what you want to do: http://www.ll.mit.edu/publications/journal/pdf/vol04_no2/4.2.5.neuralnetwork.pdf

Fortunately this is a very active field and much of the code necessary for this type of problem is readily available online.

Once you are able to identify the glasses in the images you can perform further analysis from there.


Spread spectrum receiver I want to design a receiver for a spread spectrum system. I know that two Pseudo sequence with two different length must be used in the receiver, but I can' t manage this with common receivers(such as integrate and dump and so on...) , can this be done?
Regards.

Just despread the same way with two different integrate and dump times for your two symbol types with appropriate blanking on each between longer and shorter code sequence periods.

I would be tempted to do all the acquisition and tracking on the longer symbol and then simply despread the shorter symbols using the synchronization from the longer one— with that approach envision a receiver for the just the longer sequence with 62 chip gaps between each symbol.

Alternatively if each code repeats at the same location in each symbol, given the relatively short sequences, a sliding complex correlator using I and Q FIR filters with the sequence as the coefficients (in reverse order) using two samples per chip would be a very simple approach —- the complex I and Q output of the filter would have a spike in time after every alignment with those symbols (0 or 180 degrees once aligned depending on the data symbol and given I and Q would have phase to provide for carrier recovery).  Detect the spike with a threshold detector on the magnitude for acquisition, track the complex peak for carrier tracking and track the two samples on either side of the peak for early-late timing tracking. The peak sample will give you the same performance as an integrate and dump of the multiplication of the received sequence with the output of a code generator but given your peculiar waveform structure this may be a lot simpler.

Once all timing is established, a single 31 tap FIR filter on I can simply demodulate the other two symbols.

Summarizing this FIR approach:

One 126 tap FIR filter on I  for carrier and timing recovery/tracking and demod of 63 chip symbols.

One 126 tap FIR filter on Q for carrier and timing recovery/tracking.

One 31 tap FIR filter on I for demod of the 31 chip symbols.

As @WilliamHa points out in the comments, this approach will only work if the same code sequence appears in each symbol, starting at the same location each time (beyond being inverted to designate a "1" or "0" data bit), but will not work if the code generators run continuously causing the code starting location to shift by one sample at the start of each longer sequence, and at the start of the pair of shorter sequences, due to the relationship between the two codes.  In this case I would approach it using continuously running correlators where the outputs are blanked during the other symbols, using the longer code for tracking and carrier recovery and using the synchronized clock from that to demodulate the shorter codes. Similar to the FIR filter filter approach, there would be 6 correlators for the longer code: 2 for Early (I and Q), 2 for prompt (I and Q) and 2 for late (I and Q). Each of these should be running at at least 2 samples per chip. Only one I correlator is needed for the shorter codes at one sample per chip as the carrier and time reference is established from the longer code. The added complexity will be in timing the integrate and dump outputs of the correlator such that they are blanked during the shorter code sequences, but this would be determined in the timing recover process such that accumulation starts at the estimated start of the longer code and dumps after 63 chips, and then waits 62 chips before accumulating again. This structure assumes the code generators run continuously with the same chip rate, and the correlator structure is done by multiplying the code generator output with the received signal after having been down-converted to baseband. With this approach Carrier Recovery is done by monitoring the rotation of the I and Q prompt symbols of the longer code; a carrier error is determined using a complex conjugate multiply of one symbol with the previous and then multiplying that time I prompt to strip the +1/-1 data. Clock Recovery can be done by monitoring the magnitude of the Early and Late outputs of the longer code, where clock error = Early-Late. Acquistion is done by setting a threshold on the magnitude of the longer code.

Using the cross correlation property of the FFT, all the correlators could be done with FFT's as well, including one-shot joint frequency and delay acquisition if processing allows. I explain that approach further in this post: GPS signal acquisition

prn_code= 512;%sampled at Chip2x
rx_sig % Load your signal
early = rx_sig.*prn_code(+1 circular shift);
late = rx_sig.*prn_code(+1 circular shift);
if (early-late)>1
    prn_code= prn_code(-1 circular shift )
else
    prn_code= prn_code(+1 circular shift )
end
prompt = rx_sig.* prn_code;

Am I right? @Dan Boschen


How to shift the frequency spectrum? Lets say we have a spectrum ranging from -X MHz to +X MHz. I would need to correct the frequency error in the spectrum by shifting the zero component to the middle (0 Hz).

If the output (the frequency spectrum) is calculated via FFT, as far as I know I can move the spectrum by adjusting the 'twiddle factors' (or coefficients, for complex data sine and cosine waves).

In the case of a size 1024 FFT (bin indexes from 0 to 1023), 0 Hz component should exist in bin number 511. However, due to possible frequency error the 0 hz component may actually be in bin 510 for example.

I cannot seem to find much information on this. Any help appreciated.

EDIT: Mistake in the question.

If the frequency shift that you want is a multiple of the bin spacing, as in your example, then you can easily effect the shift that you want by just rotating the FFT outputs by the number of bins that you need. In the more common case that the frequency offset is not an integer multiple of the bin spacing, then you can multiply the signal by a complex exponential function before doing the FFT.

So, if you determine that the center frequency component that you speak of is actually located at frequency $f_{offset}$ Hz in your data, and the data is sampled at rate $f_s$ Hz, then to shift the spectrum such that the component of interest is at zero frequency in the FFT output, you would do:

$$
x_{\text{shifted}}[n] = x[n] e^{-\frac{j 2 \pi f_{\text{offset}} n}{f_s}}, \ \ n = 0, 1, \ldots , N-1
$$

$$
Y_{\text{shifted}}[k] = FFT\left[x_{\text{shifted}}[n]\right]
$$

well, the simple way is , if u have used the fourier to find the spectra , and u need to know its frequency by how much is it shifted, u can do one thing ..

1) find out the impulse response from that spectra

2) convolve it with a noise

3) see the signal that u obtaine4) take its FFT just to be sure , if it matches with  the previous one

4)and see the spectra by averaging in to differnt parts... for this i can give u a algorithm in mathematica software, which is

reflect[a_] := Module[{n = Length[a]},
  RotateRight[a, Floor[n/2]]
  ]

freqAxis[len_] := Module[{},
   If[OddQ[len],
    Range[1, len] - (Ceiling[len/2.]),
    Range[1, len] - (1 + Ceiling[len/2.])
    ]
   ];

colors = {Black, Red, Blue, Brown , ColorData["Legacy", "DarkGreen"],
   ColorData["Legacy", "Goldenrod"], ColorData["Legacy", "DeepPink"],
   Cyan, Orange, Purple, ColorData["Legacy", "DeepSkyBlue"], Magenta};

specPlot[pieces_, pieceLen_, color_] :=
 Module[{data, spec, fAxis, pos},
  fAxis = freqAxis[pieceLen];
  data = Partition[Take[mysignal, pieces*pieceLen], pieceLen];
  spec = Total[Abs[Fourier[data]]^2]/pieces;
  spec = reflect[spec];
  Print["valley=", Nearest[spec, 1.0][[1]], " atPos=",
   pos = Position[spec, Nearest[spec, 1.0][[1]]][[1, 1]], " atFpos=",
   Position[fAxis, 0][[1, 1]], " atF=", fAxis[[pos]], " firstMax=",
   Max[Take[spec, Round[pieceLen/2]]], " atF=",
   fAxis[[Position[spec, Max[Take[spec, Round[pieceLen/2]]]][[1,
     1]]]], " lastMax=", Max[Take[spec, -Round[pieceLen/2]]], " atF=",
    fAxis[[Position[spec, Max[Take[spec, -Round[pieceLen/2]]]][[1,
     1]]]]];
  ListLinePlot[Transpose[{fAxis, spec}], PlotStyle -> colors[[color]],
    PlotLabel -> "N = " <> ToString[pieces], PlotRange -> All]
  ]

in this code , i have a argument to take pmsesignal , so u can use ur own signal instead of it..

i am not sure, how well i explained this, but this had worked in my case..

Cheers!


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.

I think it is just a handful way to say

$R(s,t)= E[X(s)X(t)] = \sigma^2$ if $s=t$, otherwise $0$.
If we express $R(s,t)$ as  $R(\tau)$, then if $\tau=0$, $R(\tau) = R(0) = \sigma^2$, else if $\tau \ne 0$, $R(\tau) = 0$.

The conclusion is that we can express the formula in a more compact way by using the delta function $\delta(n)$.


Why DFT-SPREAD OFDM has lower PAPR than OFDM? I just read 4G LTE/LTE-Advanced
for Mobile Broadband, and found the content as below,

  If the DFT size M equals the IDFT size N, the cascaded DFT/IDFT
  processing would obviously completely cancel each other out. However,
  if M is smaller than N and the remaining inputs to the IDFT are set to
  zero, the output of the IDFT will be a signal with “single-carrier”
  properties – that is, a signal with low power variations, and with a
  bandwidth that depends on M.

My question is that why DFT-SPREAD OFDM has lower variations in comparison with OFDM? The picture is DFT-SPREAD OFDM modulation scheme.

Thank you.

It is because the time domain symbols have lower PAPR than the frequency domain symbols.

The PAPR of OFDM is due to the 2 reasons. One is due to the IDFT(IFFT) and the other is due to the symbol type. The IDFT converts the time domain symbols to frequency domain symbols that have very high PAPR. The PAPR of 16QAM is also higher than that of QPSK. However, even when these symbol types are the same, the PAPR of DFT-SPREAD OFDM is still lower than that of regular OFDM. The reason is the (N point) IDFT effect is mostly cancelled out by the (M point) DFT. If the M=N (no zero), they are completely canceled out and the symbols would be just a time domain symbols. Even though N is much larger than M, the resulting symbol has more like time domain symbols, then the PAPR is lower. LTE up-link uses this type of OFDM that is called SC-FDMA. SC means Single Carrier because the symbols look like single carrier time domain symbols even though it has an OFDM structure. They use SC-FDMA for the up-link LTE for the lower PAPR seeking the longer battery life.

The IDFT converts the time domain symbols to frequency domain symbols that have very high PAPR
IDFT converts from frequency to time.


Selectively improve the frequency resolution of FFT (DFT) I have a digitally sampled signal(sampling frequency = 4000 Hz) which contains analog frequencies from 0 to 1000 Hz.
I wish to have a frequency resolution of 100 Hz(like 0,100,200,300 Hz .. 1000 Hz) and I have to improve the frequency resolution at the 0 to 100 Hz range by 10 Hz. That is at the frequencies 0,10,20,30,... 100.
So the frequency components I wish to detect are [0,10,20,...,100,200,300,...,1000].
In other words, I want more frequency resolution at the 0-100Hz range and less resolution at 100-1000Hz range.
How this is possible through FFT or related methods.

You can compute any single FFT result bin using a complex Goertzel  algorithm or filter.

Realize however that if you want to compute a number of bins on the order of log N or greater (such as 10 extra bins with 10 Hz bandwidth and spacing), it may be faster to do the full longer FFT and throw away all the extra FFT result bins.

So just compute 2 FFTs of 2 different lengths and ignore all the unwanted results.

You can't improve frequency resolution in a given frequency region without increasing a number of samples to be analysed. The picture is for example:

There are two input signals of the same frequency, but one is 64us duration while other is 256us duration. Their Fourier transforms are represented in the picture by blue and green curves respectively. Now we say that DFT (and also FFT) is sampling of the Fourier transform. Red line and circles in the picture are DFT bins of the short signal. They are obviously uniformly spaced, but changing coefficients of DFT matrix we can of course sample the spectrum in the interesting region more precisely. But the main idea is: the frequency resolution stays the same. It's dependent on signal duration (or number of samples under processing). It's clearly seen from the picture: green curve shows greater frequency resolution because it was calculated from more signal samples.


Converting from a phase response to a group delay response and vice-versa I'm trying to figure out how to convert any phase response to the corresponding group delay response.

Is there a way to do that, and vice-versa if possible?

The group delay of a filter is defined as minus the change in the phase response with respect to frequency.  If the phase response of a filter is $\Phi(\omega)$, the corresponding group delay $\tau_g$ is given by:

\begin{equation}
\tau_g = -\frac{d\Phi(\omega)}{d\omega}
\end{equation}

In Matlab code, the group delay of a 4th order Butterworth filter can be calculated like so:

[b a] = butter( 4, 0.25 ); % design the filter
[Hz fVec] = freqz( b, a ); % compute the frequency response and angular frequencies
phi  = unwrap(angle(Hz)); % extract phase response
taug = -diff(unwrap(angle(Hz)))/((fVec(2))); % compute group delay
figure; plot( fVec(2:end), taug );  % plot the result

The phase response must be unwrapped to avoid $\pi$ to $-\pi$ discontinuities.  Because the difference between adjacent frequency bins is the same, we need only scale the -diff(phi) by 1/fVec(2).  Finally, because the diff command produces a vector one less than its input, the group delay values are plotted against fVec(2:end) and not fVec.

Alternatively, you could use the matlab command grpdelay, though you're likely to learn more from the previous!

figure; grpdelay(b,a);

This should be:

tau_g = -diff(unwrap(angle(Hz)))./diff(fVec); % compute group delay


FFT spectrogram in log frequency space. How? I was reading a paper and the authors say that they extract FFT spectrograms and log-scale in frequency domain. I understand how the general spectrogram is computed, but how to make it log-scale in frequency domain? Or maybe there are libraries in Python, which already implement this one?

You can use either interpolation between FFT bins, or something like an MFCC triangular filter bank, except using equal delta_log spaced triangles instead of Mel-spaced triangles, or both.

With the speed of today's computers (and some optimized code) I have found a great advantage in doing Pitch Detection with a modified Logarithmic DFT instead of a FFT.

For pitch detection you are basically searching for a note's harmonics (partials) which coincidentally reside in the same frequency channels that we see on a piano's keyboard -- that is, 12 frequency channels within each octave (E, F, F#, G, G#, A, A#, B, C, C#, D, D#). The frequencies of the horizontal channels are to be logarithmically located along a vertical axis of frequency (see diagram below). So when I perform Pitch Detection (see link to my C++ source code) I construct a modified Log DFT which only probes each of the frequency channels which align to the frequencies of notes on a piano.

Since it is a DFT, we can individually adjust the center frequency for each channel so that they are logarithmically spaced and align with the frequencies of notes on a piano. And since it is a DFT, we can individually adjust the bandwidth for each of the channels by changing the number of samples used in summation -- more samples means a more narrow bandwidth for the channel. Because of the way that I gather data for my modified Logarithmic DFT, I do NOT have to apply a windowing function (nor do add and overlap) to the signal to get accurate results (number of samples in summation must be some multiple of samples in the channel's target period, usually the number samples in 17 periods).

If you look at the sonogram below (created by my C++ code), you will see a 3 second Logarithmic DFT from a guitar solo on a polyphonic mp3 recording. Notice how fine-tuned each of the horizontal channels are -- there is little bandwidth overlap, and the harmonics are clearly defined. For each note in this Log DFT we can identify their multiple harmonics extending vertically, because each harmonic will have the same time-width. Despite being a polyphonic recording, the presence of an individual note within the guitar solo is easily apparent to the naked eye.

Now being retired, I have decided to release the source code for my pitch detection engine within a free demonstration app called PitchScope Player. PitchScope Player is available on the web, and you could download the executable for Windows to see my algorithm at work on a mp3 file of your choosing. The below link to GitHub.com will lead you to my full C++ source code for Player, and you can either view and/or compile that Windows code.

To construct a Logarithmic DFT, study the DFTrowProbeCircQue class as defined in SPitchCalc.cpp

To step through a simple version of my Pitch Detection algorithm, study Make_NoteList_No_Audio() in the PitchPlayerApp.cpp file.

https://github.com/CreativeDetectors/PitchScope_Player

https://en.wikipedia.org/wiki/Transcription_(music)#Pitch_detection


STFT multiplied with Gaussian window vs STFT multiplied by rectangular window I am currently analysing my signal by looking at its spectrogram to determine the sinusoidal frequency content at the local sections of my signal overtime. I can do my spectrogram in 2 ways:

Directly get the FFT of each local sections. Same as I multiply each section by a rectangular window with a height 1 and get the FFT.
Multiply each section by a Gaussian window and get the FFT.

Apparently, I get 2 different pictures of my spectrogram. Number 2 is known as the Gabor Transform.

My question here is, what is the motivation for doing the Gabor Transform vs directly doing the FFT at each section of a signal?

Using a Gaussian window has one advantage over all other windows-- it yields the most compact time-frequency view of the spectrum. In other words, the product of the uncertainty in the time and frequency domains is minimized. This is due to the uncertainty principle of the Fourier transform.

Incidentally, in physics this property of the Fourier transform is the mathematical underpinning behind the Heisenberg uncertainty principle, which states that the position and momentum cannot simultaneously be known to infinite precision. Replace "position" with "time" and "momentum" with "frequency" in that statement and you can see how this applies to spectral analysis.

(There are other reasons not to use a Gaussian window as well, depending on what your application requires.)

Also, incidentally, Gabor was a physicist who first applied many mathematical tools from physics to signal processing problems, such as the analytic signal representation. It would not surprise me if his use of the Gaussian window in signal processing was carried over from quantum mechanics.

Here is the answer to my question. The explanations in this page is quite clear: http://www.ni.com/white-paper/4844/en/#toc2


Multiplication of two impulse function $\delta(t)\cdot \delta(t)=?$ If a impulse function is mulplied with a function $f(x)$ then the formula will be apply
$$f(x)\delta(x-a)=f(a)\delta(x-a)$$ so from this formula mulplication of two impulse function will be

  $$\delta(t)\delta(t)=\delta(0)\delta(t)$$

But,what if i take two impulse function as the approximation of standard signal (rectangular or gaussian) and mulplied them directy then i will get again an impulse function of same property.

What is wrong with this two?

The property

$$f(t)\delta(t-t_0)=f(t_0)\delta(t-t_0)\tag{1}$$

is only valid for a function $f(t)$ that is continuous at $t=t_0$. Since the Dirac delta impulse $\delta(t)$ is not a function (it is a distribution) and since $\delta(t)$ is not continuous, property $(1)$ does not hold (and does not make sense) for $f(t)=\delta(t)$.

The product $\delta(t)\cdot \delta(t)$ is undefined. Also note that the quantity $\delta(0)$ is meaningless. Since $\delta(t)$ is not a function you cannot determine its value for any value of its argument.

Well, since
$$\pi^2\delta(0)^2=2\int_0^\infty t dt-\frac1{12}$$
we have
$$\delta(0)^2=\frac2{\pi^2}\int_0^\infty t dt-\frac1{12\pi^2}$$
So, we can write $\delta(x)^2$ as
$$\delta(x)^2=\frac1{\pi^2}\int_{-\infty}^\infty |t| \cos(xt) dt-\frac{0^{|x|}}{12\pi^2}$$


Detect lanes from a filtered image I am working on lane detection system and I would like some suggestions from you as to how to detect lanes from a filtered image which I got after performing certain operations. Firstly, I should mention a few things:

The lane detection would be in a real road environment unlike here.
As of now, I am running it in real-time (O(n)) to be precise, which I would like to maintain.

The following image shows the original image and my filtered image.

As you can see, there are false positives which I want to remove and that is my main problem. Now, these are the few approaches I have thought of:

Check the intensity on the left and right hand side for desired areas. It should be high for desired area and low on both sides.
**Problem: I am finding it difficult to define 'high' and 'low' since the method has to be independent of time of the day.
One approach I thought was to check shapes but I can't think of any algorithm since those are not rectangles, those are essentially distorted parallelograms(?).

One more, I have seen people generate lines on these type of images using hough, this same question does that. If I run hough, I am not getting lines. Is the intensity of the image a problem or should I run Hough on binary mask image? I am currently doing this with MATLAB but I will be finally implementing with OpenCV to keep it real-time.

Finding shapes / lines is a good idea. Threshold tweaking is almost certainly a bad idea.

Hough transform typically expects your lines to be (1) in a binary image and (2) of one pixel width. Try to apply an edge detector (Matlab: edge) on your filtered image in order to get the boundaries of your white regions as one-pixel lines. Then Hough transform should find something.

I would also consider doing edge detection on your original image rather than your feature image. Then you can try to tweak the edge detector, rather than first running your "certain operations" and then detecting edges.

Once you have some lines you will face the problem of finding the good ones. You could look at (i.e. google for) "vanishing point detection", and read about RANSAC.

The following page might help

Road detetction

It gives results with different methods, and also explains them more or less.


Amplitude of frequency in MATLAB FFT I am trying to extract amplitude of specific frequency in Matlab FFT. Is it possible to use the abs(mag)... But I do not know in which sample to look for mag(245) should give me amplitude for the frequency of that sample...

How to extract that amplitude-magnitude for $120$ Hz using mag?

I will add simple code:

Fs = 1000;                    % Sampling frequency
T = 1/Fs;                     % Sample time
L = 1000;                     % Length of signal
t = (0:L-1)*T;                % Time vector
% Sum of a 50 Hz sinusoid and a 120 Hz sinusoid
x = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);
y = x + 2*randn(size(t));     % Sinusoids plus noise
f=(0:L-1)*Fs/L;
x=fft(y);
mag=abs(x);
mag(1)=0;
plot(f(1:L/2),(2/L*mag(1:L/2)));
title('Single-Sided Amplitude Spectrum of y(t)')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|')

The formula you are looking for is

  (desired freq / Sampling Frequency) * Length of samples = sample number

you can see this works out even by verifying the units

  (Hz / Hz )* sample = sample

depending on your exact frequencies, the result may or may not be an integer number. If it is not an integer, that means the exact frequency wasn't captured (due to Fs). For example if your result was 23.8 you would have to look at either sample 23 or sample 24 to get the data you required, or do some sort of interpolation.

Good luck

edit

as your comment specified, this doesn't work exactly, the reason is becasue matlab does 1 indexing, not 0 indexing. So our DC component 0Hz is at element 1, not element 0 but the above formulas are correct. However implementing in matlab is a little different

desired_freq = 120; %Hz
desired_index = (desired_freq / Fs) * L + 1; %+1 is because matlab uses 1 indexing not 0 indexing

again desired index may not be an integer number so you can't simply say mag(desired_index) one workaround would be mag(round(desired_index))

EDIT

I just noticed that you have already scaled your fft output.
Mag(121) will give you what you want.
If you still want the amplitude to be scaled like your figure do, scale Mag(121) too.

Just a problem with indexing.


Phase distortion of frequency warping by using Praat I am trying to do speaker conversion by frequency warping method. I modify the function Spectrum_shiftFrequencies() in Praat by scaling the FFT spectrum of a voice. When one scaling factor is multiplied to spectrum, the timbre of source voice is changed and the voice quality is good. However, using multiple scaling induces several frequency bands of echos into the synthesis audio. I guess that the multiple scaling method breaks the phase information of a spectrum.
Fig.1 single scaling and multiple scaling

Here is the code that I modify in Spectrum_shiftFrequencies() in Praat.

for (long i = 1; i <= thy nx; i++) {  //scan target frequency component
    double thyf = thy x1 + (i - 1) * thy dx;  // target freq
    double myf =  0;                 //  from source freq

    if(thyf>0 && thyf <=target_formant[0]) // 2 scaling
        myf=thyf*source_formant[0]/target_formant[0];
    else
        myf=source_formant[0] + (thyf-target_formant[0])*(my xmax-source_formant[0])/(my xmax-target_formant[0]);

        if (myf >= my xmin && myf <= my xmax) {
            double index = Sampled_xToIndex (me, myf); // index:  source freq index

            thy z[1][i] = NUM_interpolate_sinc (my z[1], my nx, index, interpolationDepth);
            thy z[2][i] = NUM_interpolate_sinc (my z[2], my nx, index, interpolationDepth);
      }

}

The real and image parts of target spectrum are interpolated from source spectrum. After checking the amplitude spectrum of synthesis audio, the contour of spectrum is changed as the expected scaling function. The phase spectrum is also changed, but that induces echo.

Does the phase spectrum require special process in this multiple scaling (nonlinear) conversion? or How to eliminate the echo ? or What is the correct phase conversion of this frequency mapping ?

Any help is appreciated. Thank you.

After research works, I find that the answer is the pitch shift. The terminology is quite confusing. The proper term should be "frequency shift" because it changes the spectrum of signal, not only pitch but also formant.

Some useful resource is listed bellow:

Pitch Shifting Using The Fourier Transform

guitar pitch shifter


Super- resolution: Warp Matrix The super resolution model is as such: Y_k = DHFX + N , where Y_k is one of the N measured low resolution image, D is the Decimation Matrix, H is the blur matrix , F is the warp matrix, N is the noise and X is the high resolution image. Lets say the X is [L *L] or we can say [L^2 * 1], F is [L^2*L^2], and H is [L^2*L^2] and D is [M_k^2*L^2] where M_k is the size of Y.
In the super resolution problem we need to find X but for that we need to find the F(warp matrix). The shifting is done to get the sub pixel information as far as I understand. I am also able to find the [3*3 ] warp matrix which transform the coordinate space but that is not required. So if anyone can help me to find the F matrix. Thank you

If H is representation of shift invariant operator then its a convolution operator. If you are implementing your algorithm in MATLAB, there is a intrinsic function with name convmtx2 which exactly do what you want.

Algorithms which warp image don't use a linear operator (a matrix) to do the job (considering large size of this operator it's not practical, even it was practical it's not an efficient method). They usually find the coordinates of warped image which generally don't match the pixel coordinates of original image then obtain pixel values using some interpolation method like nearest neighbor or bi-linear to obtain the pixel value from its neighbors.

If you want to do it by using a linear operator you have to choose an interpolation method then find the transformed coordinates and compare them to the original coordinates, then find the pixel weights regarding the interpolation method and store these weights in their corresponding location in the matrix (in nearest neighbor method give weight 1 for the nearest pixel and zero for others, in bi-linear method 4 non-zero weights for 4 pixel which are around the transformed coordinate and zero for other pixels).

Good luck.

Hey mohammad  its not convolution , with the 3*3 matrix we multiply it to the original coordinate axes to get the new coordinate axes.


Remove a known wav file from recorded file I looked at many similar questions but still cannot find an answer to my problem.

I'm playing an audio file N, while recording a man talking (b).
The problem is that the microphone gets the played sound N together with the man talking, b. so actually I have y=N*+b in my recording, where N* is the played sound with a short time delay and a bit noise (due to the recording and the phone speaker).

My goal is to totally remove the N* from the recorded file. It seems like an easy problem cause I have the original played file N, but I find it more complicated.

What I tried:

finding the delay between the N and y, invert N and add to y. Seems like it's not working and the result is even worse probably because N* has different amplitudes and a phase shift.
Tried the same thing but subtracting in the frequency domain between y and N. Still no good.

Any suggestions?

Assuming that there are no non-linear effects between the sound N source and the microphone (such as AGC), you might try to estimate the impulse response of the channel between the source for sound N and the microphone data.  The better estimate you can make of that impulse response, the more of sound N that you can remove.  It’s unlikely that the negative of the impulse response is just the data inverted. That's due to the fact that the amplifier, speaker, microphone, and room acoustics (etc.) probably have frequency dependent phase delays and/or multiple reverb/echoes.

Channel estimation in noise is a non-trivial problem.

So I actually just spent the last year working on a website that does just this, www.tuneout.io. I just launched a free Beta version. Curious to hear if it works for you.

Since SO doesn't like here's-a-website-that-exactly-solves-your-problem answers, I'll address some points brought up by other commenters. The problem is mainly the channel distortion on the audio (room acoustics, speaker distortion, etc.). Phase inversion only really works if the sound is very clean, with little to no reverberation, the speaker and mic have flat frequency responses, and the volume is constant. Doing the subtraction in the frequency domain is actually mathematically the same as subtraction in the time domain (assuming no windowing). So it's expected you'll get the same result. Another issue is that the song in the recording will often be a different length than your audio file of the song. This is because real-life devices all have slightly different playback speeds due to slight hardware differences. So, even if you find the exact delay using finddelay(), if there is "time drift" in the file, the delay will also drift, quickly rendering phase inversion useless. Add to this you could be starting with a file of slightly different length due to compression (mp3s can be several seconds longer than a wav file sometimes).

My algorithm uses a special blend of spectral subtraction and machine learning to address all of these problems in a robust (hopefully) fully automated way. Hope it works for you!


What is the physical significance of negative frequencies? This has been one of the holes in my cheddar cheese block of understanding DSP, so what is the physical interpretation of having a negative frequency?

If you have a physical tone at some frequency and it is DFT'd, you get a result in both the positive and negative frequencies - why and how does this occur? What does it mean?

Edit: Oct 18th 2011. I have provided my own answer, but expanded the question to include the roots of why negative frequencies MUST exist.

Negative frequency doesn't make much sense for sinusoids, but the Fourier transform doesn't break up a signal into sinusoids, it breaks it up into complex exponentials (also called "complex sinusoids" or "cisoids"):

$$F(\omega) = \int_{-\infty}^{\infty} f(t) \color{Red}{e^{- j\omega t}}\,dt$$

These are actually spirals, spinning around in the complex plane:

(Source: Richard Lyons)

Spirals can be either left-handed or right-handed (rotating clockwise or counterclockwise), which is where the concept of negative frequency comes from.  You can also think of it as the phase angle going forward or backward in time.

In the case of real signals, there are always two equal-amplitude complex exponentials, rotating in opposite directions, so that their real parts combine and imaginary parts cancel out, leaving only a real sinusoid as the result.  This is why the spectrum of a sine wave always has 2 spikes, one positive frequency and one negative.  Depending on the phase of the two spirals, they could cancel out, leaving a purely real sine wave, or a real cosine wave, or a purely imaginary sine wave, etc.

The negative and positive frequency components are both necessary to produce the real signal, but if you already know that it's a real signal, the other side of the spectrum doesn't provide any extra information, so it's often hand-waved and ignored.  For the general case of complex signals, you need to know both sides of the frequency spectrum.

This has turned out to be quite the hot topic.

After reading the rich multitude of good and diverse opinions and interpretations and
letting the issue simmer in my head for sometime, I believe I have a physical interpretation of the phenomenon of negative frequencies. And I believe the key interpretation here is that fourier is blind to time. Expaning on this further:

There has been a lot of talk about the 'direction' of the frequency, and thus how it can be
+ve or -ve. While the overarching insights of the authors saying this is not lost, this
statement is nontheless inconsistent with the definition of temporal frequency, so first we
must define our terms very carefully. For example:

Distance is a scalar (can only ever be +ve), while displacement is a vector. (ie, has
direction, can be +ve or -ve to illustrate heading).
Speed is a scalar (can only be +ve), while velocity is a vector. (ie, again, has
direction, and can be +ve or -ve).

Thus by the same tokens,

Temporal Frequency is a scalar, (can only be +ve)! Frequency is defined as number of
cycles per unit time. If this is the accepted definition, we cannot simply claim that it is going in 'a different direction'. Its a scalar after-all. Instead, we must define a new term - the vector equivalent of frequency. Perhaps 'angular frequency' would be the right terminology here, and indeed, that is precisely what a digital frequency measures.

Now all the sudden we are in the business of measuring number of rotations per unit time,
(a vector quantity that can have direction), VS just the number of repititions of some physical oscillation.

Thus when we are asking about the physical interpretation of negative frequencies, we are
also implicitly asking about how the scalar and very real measures of number of
oscillations per unit time of some physical phenomenon like waves on a beach, sinusoidal AC
current over a wire, map to this angular-frequency that now all the sudden happens to have direction, either clockwise or counterclockwise.

From here, to arrive at a physical interpretation of negative frequencies two facts need to
be heeded. The first one is that as Fourier pointed out, an oscillatory real tone with
scalar temporal frequency, f, can be constructed by adding two oscillatory complex tones,
with vector angular frequencies, +w and -w together.

$$
\cos(\omega_0 t)=\frac{e^{\jmath \omega_0t}+e^{-\jmath \omega_0 t}}{2}
$$

Thats great, but so what? Well, the complex tones are rotating in directions opposite to
each other. (See also Sebastian's comment). But what is the significance of the 'directions' here that give our angular frequencies their vector status? What physical quantity is being reflected in the direction of rotation? The answer is time. In the first complex tone, time is travelling in the +ve direction, and in the second complex tone, time is travelling in the -ve direction. Time is going backwards.

Keeping this in mind and taking a quick diversion to recall that temporal frequency is the
first derivative of phase with respect to time, (simply the change of phase over time), everything begins to fall into place:

The physical interpretation of negative frequencies is as follows:

My first realization was that fourier is time-agnostic. That is, if you think about it,
there is nothing in fourier analysis or the transform itself that can tell you what the
'direction' of time is. Now, imagine a physically oscillating system (ie a real sinusoid
from say, a current over a wire) that is oscillating at some scalar temporal-frequency, f.

Imagine 'looking' down this wave, in the forwards direction of time as it progresses. Now
imagine calculating its difference in phase at every point in time you progress further.
This will give you your scalar temporal frequency, and your frquency is positive. So far so
good.

But wait a minute - if fourier is blind to time, then why should it only consider your wave
in the 'forward' time direction? There is nothing special about that direction in time. Thus by symmetry, the other direction of time must also be considered. Thus now imagine 'looking' up at the same wave, (ie, backwards in time), and also performing the same delta-phase calculation. Since time is going backwards now, and your frequency is change-of-phase/(negative time), your frequency will now be negative!

What Fourier is really saying, is that this signal has energy if played forward in time at
frequency bin f, but ALSO has energy if played backwards in time albeit at frequency bin
-f. In a sense it MUST say this because fourier has no way of 'knowing' what the 'true'
direction of time is!

So how does fourier capture this? Well, in order to show the direction of time, a rotation of some sort must be employed such that a clockwise roation dealinates 'looking' at the signal in the forward arrow of time, and a counterclockwise roation dealinates 'looking' at the signal as if time was going backwards. The scalar temporal frequency we are all familiar with should now be equal to the (scaled) absolute value of our vector angular frequency. But how can a point signifying the displacement of a sinusoid wave arrive at its starting point after one cycle yet simultaneously rotate around a circle and maintain a manifestation of the temporal frequency it signifies? Only if the major axes of that circle are composed of measuring displacement of this point relative to the original sinusoid, and a sinusoid off by 90 degrees. (This is exactly how fourier gets his sine and cosine bases the you project against every time you perform a DFT!). And finally, how do we keep those axes seperate? The 'j' guarantees that the magnitude on each axis is always independant of the magnitude on the other, since real and imaginary numbers cannot be added to yield a new number in either domain. (But this is just a side note).

Thus in summary:

The fourier transform is time-agnostic. It cannot tell the direction of time. This is at the heart of negative frequencies. Since frequency = phase-change/time, anytime you take the DFT of a signal, fourier is saying that if time was going forwards, your energy is located on the +ve frequency axis, but if your time was going backwards, your energy is located on the -ve frequency axis.

As our universe has shown before, it is precisely because Fourier does not know the direction of time, that both sides of the DFT must be symmetric, and why the existence of negative frequencies are necessary and in fact very real indeed.


Choosing the right overlap for a window function On Wikipedia there is quite a list of window functions.  However, when you are doing a lot of STFTs, you often have to have some overlap to get better results.

Generally an overlap factor of 2 (i.e., 50% overlap) works fine (I'm doing analysis & resynthesis), but for some windows it isn't good (I currently have this experience with the Blackman window).

I tried my best but I couldn't find a resource that would list the "good" overlap factors for common and less common windows.

Please refer to this paper on the optimum overlap percentage for the Blackman-Harris window, which is derived to be 66.1%. It has a lot of useful information on spectral analysis and windows.

Maybe look for the keywords ,,shift overlap add property'' or ,,constant overlap add property'' - Some window functions have the advantage to sum up to a constant function if taken in the proper distance, as was remarked here.


Correlated fading channel: transmit and receive antenna correlation matrices I am trying to create a non-i.i.d channel: correlated Rayleigh fading channel.
$$ H= R_H \cdot H_o,$$
where $H_o$ is i.i.d. channel ( uncorrelated) and $ R_H$ is a channel correlation matrix, is defined as
$$ R_H= R_T^T \otimes R_R,$$
with $R_R$ receive antenna and  $ R_T$ transmit antenna correlation matrices.

My question is how to find $R_T$ and $R_R$, if numbers transmit and receive antennas are given only?
Can I assume them as identity matrices?

$R_H$ could be anything other than identity and you would get correlated fading with the definitions you have in place in the question. Defining a kronecker matrix product is one way.

$R_H$ should be $N_R \times N_R$ matrix, you can make the kronecker product to result in a non identity matrix of this size.


I often see here formula expressed in term of $z$. But what is $z$? While searching resources for generating pink noise (and with your help in the comments and answers of other questions), I came to such kind of formula:

  $$ H(z) = { .041 - .096z^{-1} + .051z^{-2} - .004z^{-3} \over{}  -2.495z^{-1} + 2.017z^{-2} - .522z^{-3}} $$

  (from https://dsp.stackexchange.com/a/376/46389)

What is $z$?

I know it should be obvious for anyone familiar with digital signal processing but for people coming from a different field, it is not obvious to know what this represents and how it relates to the signal amplitude or frequency or any other "tangible" parameter.

Edit: amusingly enough, the site suggested me the z-tranform tag for this question. I suspect I should consider that as a clue.

$z$ is interpreted as the time advance operator. $z^{-1}$ is the time delay operator.

So for a difference equation like
$$ y[n]=a x[n]+b x[n-1] $$
in the $z$ domain
$$ Y(z)=(a+ b z^{-1}) X(z)$$

$z = e^{j\omega}$ is a convenient substitution for a complex valued function.


What mathematical function could be used to simulate this sound effect? I have this sound effect which by the way if anyone could give me its name, I would really appreciate it. And after loading it into matlab/octave (don't worry this isn't a programming question) using the following commands:
[x, fs] = audioread("Path To File");
plot(length(x(:,1)),x)
sound(x,fs) % To hear the sound

I get the following:

Could anyone tell me what mathematical function I could use to simulate this?

I guess you want to synthesize that sound, i.e. create a synthetic signal to be as close as possible to the original sound. Potentially, you can do this by creating a superposition of sine waves with varying amplitudes and phases. So you have to know the parameters for those sine functions.
From what I can see and hear, the sound has two signals, which are very similar (both visually and to ear). Moreover, if you look at the frequency spectrum, you will see, that they are quite similar as well. However, the spectrum seems quite populated with different frequencies, which will make it difficult to reproduce as a synthetic function. At least if using only the provided original sound file.
Here is the spectrum plot:

And here is a python code to make those observations. (The file with sound is downloaded to the folder of the python script.)
import matplotlib.pyplot as plt
import soundfile as sf
import sounddevice as sd

# get file
filename = "SoundEffect.wav"
x, f = sf.read(filename)
x1 = x[:,0]     # the first signal
x2 = x[:,1]     # the second signal

# play and show the original sound in time domain
sd.play(x, f)
status = sd.wait()
plt.plot(x)
plt.show()

# play the first signal
sd.play(x1, f)
status = sd.wait()
# play the second signal
sd.play(x2, f)
status = sd.wait()

# show the signals in frequency domain
plt.magnitude_spectrum(x1, Fs = f)
plt.magnitude_spectrum(x2, Fs = f)
plt.show()

Bear in mind, I am not familiar with how sound effects are done. Most probably, this sound effect is a more simple superposition of certain signal shapes. However, unless they are already known, it is hard to "reverse-engineer" them from the sampled version of the signal.


How to Extract Nonorthogonal PCA Principal Components My problem is essentially a 'blind source separation' problem. I have 3 non-orthogonal sources (or basis functions) and $N$ random linear combinations (mixes) of said sources. My problem is to obtain the sources from the mixes.

Figure A shows the sources, B shows the mixes.

Approaches taken:

$\textbf{PCA}$ - I tried PCA on the mixes (fig.C), but the issue is that PCA will only give orthogonal bases, while my sources/bases are non-orthogonal. This issue with PCA is shown in figure D, where the data is clearly described by 2 non-orthogonal basis, but PCA (the solid lines) cant reconstruct them!
$\textbf{Factor Rotation}$ - I tried applying some solutions form Factor Analysis (not my forte). The promax rotation (matlab: nw = rotatefactors(cov,'method','promax') ) is shown in fig. D with the dashed lines. As far as I can tell, factor rotation works with the principal components, not the original matrix and thus I have no idea how it could reconstruct the right basis. I think this only works with factor matrices, not with generic ones...
$\textbf{ICA}$ - I tried overdetermined ICA with the fastICA algorithm (also not my forte, but I think I understand it better). I was hoping this would work since independent components (ICs) are non orthogonal. The solution is shown in fig.E. While the ICs are indeed nonorthogonal, they are NOT my original sources :-(

Any other potential tips or leads or solutions would be greatly appreciated.

More mathematical Detail

The 3 source signals are themselves linear combos of the 5 laguerre basis functions. Thus, the problem can be formulated as such:

I take l=5 laguerre basis functions of length m=30 to get B, an [l,m] matrix
I randomly mix then is my x=3 sources, S: S=B*M1, where M1 is a random [l,x]=[5,3] mixing matrix. S is displayed in fig.A
I mix S n=100 times to get my mixed signals, K (fig.B), where K=S*M2, a [x,n]=[3,100] randomly generated mixing matrix.
I have the original basis, L, & the mixed matrix, K. K=S*M2=B*M1*M2. My goal is to recover S.

NOTE: I understand that this may seem impossible since there are many possible S's to choose from (ie the factorizatio S*M2=B*M1*M2 is not unique), however, the true S, will dramatically decrease model complexity by having only x=3 rather than l=5 basis functions which can exactly reproduce all the kernels.

Thanks!

Let's say x,y and z are your basis vectors and you mix them to get what you are calling your mixed sources. If your basis vectors were (x-y), y, and (z+y) then you could mix them in a different way and end up with the same mixed sources.

Basically, there are many ways to come up with a set of basis vectors for your signals. You need to say what kind of basis you want. It is not enough to say that the basis vectors are non-orthogonal.

A map from a weights vector $w$ to a ,,mixed signal'' $S$ would be just your 3 prototype signals as a matrix

$A=\left[\matrix{x & y & z}\right]$

so that a linearly combined signal could be written as

$S=Aw$

A ,,pseudoinverse'' to this matrix $A$ (in a few words: a matrix $A'$ so that $A'A=\mathbf{1}$) would map that signal back to $w$:

$A'S=w$

The problem is then finding such a pseudo-inverse. However, this problem is well-known and implementations are available.


Difference between discrete time fourier transform and discrete fourier transform I have read many articles about DTFT and DFT but am not able to discern the difference between the two except for a few visible things like DTFT goes till infinity while DFT is only till N-1. Can anyone please explain the difference and when to use what? Wiki says

  The DFT differs from the discrete-time Fourier transform (DTFT) in
  that its input and output sequences are both finite; it is therefore
  said to be the Fourier analysis of finite-domain (or periodic)
  discrete-time functions.

Is it the only difference?

Edit:
This article nicely explains the difference

alright, i'm gonna answer this with an argument that "opponents" to my rigid nazi-like position regarding the DFT have.
first of all, my rigid, nazi-like position:  the DFT and Discrete Fourier Series is one-and-the-same.  the DFT maps one infinite and periodic sequence, $x[n]$ with period $N$ in the "time" domain to another infinite and periodic sequence, $X[k]$, again with period $N$, in the "frequency" domain.  and the iDFT maps it back.  and they're "bijective" or "invertible" or "one-to-one".
DFT:
$$ X[k] = \sum\limits_{n=0}^{N-1} x[n] e^{-j 2 \pi nk/N} $$
iDFT:
$$ x[n] = \frac{1}{N} \sum\limits_{k=0}^{N-1} X[k] e^{j 2 \pi nk/N} $$
that is most fundamentally what the DFT is.  it is inherently a periodic or circular thing.
but the periodicity deniers like to say this about the DFT.  it is true, it just doesn't change any of the above.
so, suppose you had a finite-length sequence $x[n]$ of length $N$ and, instead of periodically extending it (which is what the DFT inherently does), you append this finite-length sequence with zeros infinitely on both left and right.  so
$$ \hat{x}[n] \triangleq \begin{cases}
x[n] \qquad & \text{for } 0 \le n \le N-1 \\
\\
0 & \text{otherwise}
\end{cases} $$
now, this non-repeating infinite sequence does have a DTFT:
DTFT:
$$ \hat{X}\left(e^{j\omega}\right) = \sum\limits_{n=-\infty}^{+\infty} \hat{x}[n] e^{-j \omega n} $$
$\hat{X}\left(e^{j\omega}\right)$ is the Z-transform of $\hat{x}[n]$ evaluated on the unit circle $z=e^{j\omega}$ for infinitely many real values of $\omega$.
now, if you were to sample that DTFT $\hat{X}\left(e^{j\omega}\right)$ at $N$ equally spaced points on the unit circle, with one point at $z=e^{j\omega}=1$, you would get
$$ \begin{align}
\hat{X}\left(e^{j\omega}\right)\Bigg|_{\omega = 2 \pi\frac{k}{N}} & = \sum\limits_{n=-\infty}^{+\infty} \hat{x}[n] e^{-j \omega n} \Bigg|_{\omega = 2 \pi\frac{k}{N}} \\
& = \sum\limits_{n=-\infty}^{+\infty} \hat{x}[n] e^{-j 2 \pi k n/N}  \\
& = \sum\limits_{n=0}^{N-1} \hat{x}[n] e^{-j 2 \pi k n/N}  \\
& = \sum\limits_{n=0}^{N-1} x[n] e^{-j 2 \pi k n/N}  \\
& = X[k]  \\
 \end{align} $$
that is precisely how the DFT and DTFT are related.  sampling the DTFT at uniform intervals in the "frequency" domain causes, in the "time" domain, the original sequence $\hat{x}[n]$ to be repeated and shifted by all multiples of $N$ and overlap-added.  that's what uniform sampling in one domain causes in the other domain.  but, since $\hat{x}[n]$ is hypothesized to be $0$ outside of the interval $0 \le n \le N-1$, that overlap-adding does nothing.  it just periodically extends the non-zero part of $\hat{x}[n]$, our original finite-length sequence, $x[n]$.

DFT:
$$
X[k]=\sum_{n=0}^{N−1} x[n] e^{−j2\pi nk/N}
$$
ITS INVERSE WILL BE :
$$
x[n]=\frac{1}{N} \sum_{k=0}^{N−1} X[k] e^{j2\pi nk/N}
$$


Check the stability of the system, $y[n]=u[n]$ How do I check the following system
$$
y[n]=u[n]
$$
is BIBO stable or not ?

$u[n]$ is the unit step function

My Attempt:

For the BIBO stability, the necessary and sufficient condition is $$\sum_{n=-\infty}^{+\infty} |h[n]|<\infty$$.

Let $x[n]=\delta[n]$, so $y[n]=h[n]=u[n]$
$$
\sum_{n=-\infty}^{+\infty} |h[n]|<\infty=\sum_{n=0}^{+\infty} |u[n]|=1+1+1+......=\infty
$$

which proves it is an unstable system. Is it the right way to approach the problem or do I need to first write $u[n]$ in terms of the delta function ?

(Long Comment) Sorry guys, I did not get why that system created so many conceptual problems to you. For instance, consider the unforced system $y(n)=0$. If I define $\pi:\{0\}\times\mathbb Z\to\{0\}$, $\pi(y,n)=0$ then $(\{0\},\mathbb Z,\pi)$ defines a dynamical system on the state space $\{0\}$. In fact $\pi$ is continuous, and it satisfies the identity axiom ($\pi(x,0)=x$ for all $x\in\{0\}$) and the group axiom ($\pi(\pi(x,n_1),n_2) = \pi(x,n_1+n_2)$ for all $(x,n_1),(x,n_2)\in\{0\}\times\mathbb Z$).

I agree that this is a degenerate case. I guess I can extend the state space to $\mathbb R$ by defining the system in a different way. For instance, let
$$
x(n+1)=0\qquad (1)
$$
be a ordinary difference equation and define $\pi$ as the flow of (1). Then $(\mathbb R,\mathbb N,\pi)$ defines a dynamical system on $\mathbb R$. Take the forced system
\begin{align*}
x(n+1) &= 0\cdot x(n) + 0\cdot u(n)\\
y(n) &= 0\cdot x(n) + 1\cdot u(n)\qquad\qquad (2)
\end{align*}
That makes sense to me. It is a LTI system with zero relative degree. To check BIBO stability let take a reference input $\bar u:\mathbb N\to\mathbb R$ and a perturbation $\delta u:\mathbb N\to\mathbb R$. Let $u(n)=\bar u(n)+\delta u(n)$.
Let $\bar y:\mathbb N\to\mathbb R$ be the solution corresponding to the input $\bar u$ and let $y:\mathbb N\to\mathbb R$ be the one corresponding to the input $u$. Then
\begin{align*}
\|y(n)-\bar y(n)\| = \|u(n)-\bar u(n)\|=\|\delta u(n)\|
\end{align*}
So trivially, for any $\epsilon>0$ there exists $\delta_\epsilon>0$ such that
$\|\delta u(n)\|\le \delta_\epsilon\implies \|y(n)-\bar y(n)\|\le\epsilon$. You can see this also as a continuity property of the solutions. In fact, you can define the space $Y$ of all the bounded output trajectories $y:\mathbb N\to\mathbb R$ and the space $U$ of all the bounded input ones. Let endow $Y$ with the norm $\|y\|:=\sup_{n\in\mathbb N}|y(n)|$ and $U$ with the norm $\|u\|:=\sup_{n\in\mathbb N}|u(n)|$. Then the system (2) definest an operator $T:U\to Y$, and asking for BIBO stability is the same as asking for the continuity of $T$ in the sup norms.


Getting frequencies out of a FFT I use this FFT and wrote a short program to test it.

package test_FFT;
import org.apache.commons.math3.transform.FastFourierTransformer;
import org.apache.commons.math3.transform.DftNormalization;
import org.apache.commons.math3.transform.TransformType;

public class FFT2 {

public static void main(String[] args) {
    double[][] array=new double[2][8];
    //An array of complex numbers with real part array[0][*]
    //and complex part array[1][*]
    //The output of the transformation will be saved in the input array
    array[0][0]=5.0;
    array[0][1]=2.0;
    array[0][2]=3.0;
    array[0][3]=4.0;
    array[0][4]=5.0;
    array[0][5]=6.0;
    array[0][6]=7.0;
    array[0][7]=9.0;
    FastFourierTransformer.transformInPlace(array,DftNormalization.STANDARD,TransformType.FORWARD);

    for(int i=0;i<8;i++){
        System.out.println("real "+array[0][i]);
        System.out.println("complex "+array[1][i]);
    }

}

}

In the output array I get amplitudes of sin and cos functions. The information about the frequencies should depend on the position within the output array.
After some research on this page, I still don't understand how to calculate frequencies out of array positions.

I learned that there are many flavours of how to perform a FFT. Has anyone of you detailed knowledge on how calculate frequencies out of output of the FFT I use? A code sample computing frequencies for the output of the example below would be greatly appreciated. Thank you

At this moment i'm not hos i can't check my answer...
You say

  In the output array I get amplitudes of sin and cos functions.

Well that is you frequency array. Probably not the way want it. I suspect youre looking for the magnitude and maybe the phase.

the magnitude is quite easy.. it is the $abs(array (real + imag)$.
the phase $ \phi = \arctan \frac{imag}{real} $

in this way it is common to show the half the frequency array because of it's symmetry.
if you plot $array[0]$ to $array[N/2]$ you shoult get your frequecy array from $ \frac{1}{N}$ to $\frac{F_s}{2}$