How can I convert audio to a series of notes using pitch detection? I have a good pitch detection system set up, and I would like to return a series of notes given an array of audio samples.

My current approach is as follows: I have a moving window across the audio signal, and I calculate the pitch of each window. Afterwards, I segment the audio into different notes by detecting the silent regions (i.e. where the pitch detector returns null). I then simply take the average of each note region.

Unfortunately, this has not been giving me such good results. The pitch detector does seem pretty accurate, but the issue is it doesn't seem to segment the notes very well. It only really works when I leave a long pause between each note when I record the audio. I would like some way for it to detect a change in notes without having to rely on a large silent region.

Any ideas would be greatly appreciated!

This should work using a threshold on the by-bin-difference of the magnitude spectrum.

Calculate the abs fft for current window and normalize it.
Do a bin-wise difference with the normalized abs fft of the last window.
Sum up the bin-wise difference and compare it to a threshold.

Explanation: A change of note means a change of energy distribution over the fft bins. So the sum of bin-wise difference should be high, if a change of note has occured. Normalisation is to suppress false detection of note change due to change of volume.


Clarification Needed: Power Spectral Density of Real Signal For a real signal, the power spectrum will have Hermitian symmetry. If I consider the one-sided spectrum (as negative frequencies, or anything above the Nyquist frequency are redundant) must I then multiply my power spectrum by a factor of 2? I have seen this done but I can't understand why the "redundant", symmetrical portion of the spectrum should factor into the power density? Maybe this is a matter of convention? Please correct my assumptions if necessary.

This stems from Parseval's theorem. For any signal, $\int_{-\infty}^{+\infty}|x(t)|^{2}\mathrm dt = \int_{-\infty}^{+\infty}|X(f)|^{2}\mathrm df$ and hence for a real signal, the one sided PSD must be multiplied by a factor of 2 because of Symmetry. Intuitively, the Fourier Analysis of any signal provides us with a spectrum which has Negative frequencies too because of the fact that Sines and Cosines can be constructed using Euler's theorem of using Complex Exponentials of positive and negative frequencies; and it so happens for any real signal that all you need are complex exponentials which are conjugates of each other to construct them.

The negative portion is not exactly "redundant". The signal still owns those "negative" frequencies and hence anything present there accounts to its energy.


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.


How to Generate the Following Image Effect (Modulation)? I was given the original image and I need to produce Image 1. Does anybody know how I would go about this? I am on MATLAB, but I am more interested in the theory behind this.

Original Image

**Processed Image **

NOTE How would you achieve this through convolution of the ffts of the image and a sine wave along the vertical axis?

The image seems to be modulated, namely multiplied by Harmonic Signal along its columns.

This is the MATLAB Code:

mInputImage = im2double(imread('2vE0G.png'));
numRows = size(mInputImage, 1);
numCols = size(mInputImage, 2);

signalFreq = floor(numRows / 1.25);

vModulationSignal = 1 + (0.5 * cos(2 * pi * (signalFreq / numRows) * [0:(numRows - 1)].'));

mOutputImage = bsxfun(@times, mInputImage, vModulationSignal);

figure();
imshow(mOutputImage);

This is the result:

You can play with the parameters.

Enjoy...


Why must the angle contain $\pi$ for $\cos$ be periodic? The $\cos(\frac{n}{6})$ is aperiodic, and $\cos(\frac{n\pi}{6})$ is periodic, why? What role does the $\pi$ play? By the way, $n$ must be integer

I have seen this question,Why is cos(n/6) aperiodic? but the answer just show that it is not of the form $2\pi(mN)$,but it did't explain why should that $\pi$ exist,and the comment below said that we can't find the next larger (integer) value of n such that $\cos(\frac{n}{6})=1$ after $n=0$, however, yes we can, the $n=360 \times 6$ can let $\cos$ value back to the $0$ again. So, I want to ask why must the $\pi$ exist to let $\cos$ be a periodic?

I will reuse a former answer of mine too: Proof of complex conjugate symmetry property of DFT. It relates to how sines and cosines can be defined. One answer is: from the exponential, and thus derives $\pi$. There exist other constructions, this one is (imho) elegant.

The periodicity of the cosine comes from the fact it is defined as the real part of the cisoid, or complex exponential, $e^{ix}$.
In W. Rudin's Real and Complex analysis (very first pages, 1 to 3 of the prologue), $e$ comes first, and $\pi$ appears subsequently.

  One first defines for any complex $z$:

  $$e^z=\sum_{n=0}^{+\infty} \frac{z^n}{n!}$$

  which is an absolutely convergent series. It is its own derivative. And then you get some other results, like
  $e^z$ is never equal to zero. But the two most striking ones are:

    There exists a positive number $\pi$ such that $e^{\pi j/2} = j$ and such that $e^z = 1$ if and only if $z/(2\pi j )$ is an integer.

    $e^z$ is $2\pi j$ periodic.

From this, you define the sine and the cosine as the imaginary and real parts. The proof is quite interesting. It for instance remarks that $\cos 2 <0$, and since  $\cos 0 =1$, by continuity, there should exists a constant at which the cosine vanishes:

This should answer to Why does $\pi$ exist, with its role unveiled, and why the cosine is periodic.

The answer to the "next integer" is given by Matt's comment.


Matched filter and cross correlation in digital processing I am very new to signal processing, and a bit confused with the matched filter.
Assume that I have a time series and a specific waveform I need to identify in it. By definition:

  a matched filter is obtained by correlating a known signal, or
  template, with an unknown signal to detect the presence of the
  template in the unknown signal.

So I take the window from the sequence, generate the test signal with the same sampling rate, cross-correlate them, find the maximum and threshold it.
Is cross-correlating my measured time sequence with the generated waveform a matched filter (for example, using xcorr in matlab), or there is a specific procedure I must follow to produce a "more correct" matched filter.
Thank you.

I'll illustrate with an example: detecting a rectangular pulse in noise, in Matlab.

Let's start with defining the signal that contains the pulse.

fs = 1/1000;
t = 0:fs:1;
s = [zeros(1,400), ones(1,100), zeros(1,501)];

First let's do detection in noise with small power. We calculate the received (noisy) signal:

np = 0.1;  % noise power (variance)
n = sqrt(np)*randn(1,1001); % noise signal
r = s + n; % received signal
plot(t,r); title('Received signal, small noise power');

Now let's run the correlation:

c = xcorr(r,s);
plot(c); title('Correlation, small noise power');

You can see how the correlation produced a large output, indicating that the pulse was there. What is amazing is that the correlator can find the pulse even when it is impossible for one to see it. Let us increase the noise power and repeat the process.

np = 4;
n = sqrt(np)*randn(1,1001);
r1 = s + n;  % received signal with large noise

As you can see, the correlator has an incredible ability to detect signals that are embedded in strong noise.

One thing you can do for yourself is run the correlation when the pulse is not present. The signal will be hard to distinguish from the noisy pulse shown above, but you'll see that the correlation does not produce a sharp maximum.

So, the implementation is quite easy; what is not so easy to do is:

Prove that this technique produces the best possible signal to noise ratio at the correlation's peak.
Determine the optimum threshold given some desired constraints on the probability of wrong detection (false positive or false negative).


Hamming window in frequency domain and back I can't figure out the scaling and slight artifact in fourth plot.

N             = 32;
window        = hamming(N);

A             = fft(window, 256);

A_ifft        = (N/2) * ifft(A)

subplot(411)
plot(window);
ylabel('w(n)')
title('Hamming window')

magnitude     = abs(fftshift(A));
subplot(412)
plot(magnitude);
ylabel('W(j\omega)')

response      = 20*log10(magnitude);
subplot(413)
plot(response);
ylabel('W(j\omega) dB')

subplot(414)
plot(abs(fftshift(A_ifft)));
ylabel('IFFT(W(j\omega))')

Your initial plot doesn’t show the full input to the FFT, which is 256 in length, and IIRC zero-padded.


How to apply band-pass filter and get the filtered back from frequency to time domain in python? I know this is basic for signal processing, but, I am not sure what is wrong about my approach. I have a signal that behaves as damped sine signal with a sampling frequency of 5076Hz and 15,000 number of samples. I found from the following website how to convert a signal from a time domain to frequency domain and managed to get the FFT and frequency values. The code can be found below the link:
Machine Learning with Signal Processing Techniques
def get_fft_values(y_values, T_s, N, f_s):
    f_values = np.linspace(0.0, 1.0/(2.0*T), N//2)
    fft_values_ = np.fft.rfft(y_values)
    fft_values = 2.0/N * np.abs(fft_values_[0:N//2])
    return f_values, fft_values

I managed to get the frequency and FFT values. However, I need to implement filters to remove some noise from the signal, so, I created the following functions to implement the filter part:
def butter_bandpass(lowcut, highcut, fs, order):
    nyq = 0.5 * fs
    low = lowcut / nyq
    high = highcut / nyq
    b, a = butter(order, [low, high], btype='bandpass', output='ba')
    return b, a

def butter_bandpass_filter(data, lowcut, highcut, fs, order):
    b, a = butter_bandpass(lowcut, highcut, fs, order=order)
    y = filtfilt(b=b, a=a, x=data)
    # y = lfilter(b=b, a=a, x=data)
    return y

I know that I would need to implement the following steps:

Convert to the frequency domain
apply a bandpass filter to get rid of frequencies you don't care
about
convert back to the time domain by inverse Fourier transform

So, I created the following inverse transform function, but, I can't get the filtered signal back and the amplitudes don't almost match the original signal.
(For my case, I need to resample)
def get_ifft_values(fft_values, T, N, f_s):
    # Time axis:
    N = 9903
    S_T = 1 / S_F
    t_n = S_T * N  # seconds of sampling
    # Obtaining data in order to plot the graph:
    x_time = np.linspace(0, t_n, N)
    ifft_val = np.fft.irfft(fft_values, n=N)
    y_s, x_time = scipy.signal.resample(x=ifft_val, num=N, t=x_time)
    return x_time, y_s

My Approach and here are the results from the signal:
##### Converting the signal into fft:
f_val, fft_val = get_fft_values(y_values=y, T=S_T, N=N, f_s=S_F)
# Applying bandpass filter:
fft_filt_val = butter_bandpass_filter(data=fft_val, lowcut=50, highcut=600, fs=S_F, order=2)
# Applying the inverse transform of the frequency domain:
x_time, y = get_ifft_values(fft_values=fft_filt_val, T=S_T, N=N, f_s=S_F)

FFT of the original signal:

Filtered FFT of the original signal:

Converted Signal from Filtered FFT:

Without Applying the bandpass filter:

What am I doing wrong here?

You shouldn't perform the filtering operation on the FFT data, but on the original time domain data. The functions filtfilt and lfilter both take time domain data as their inputs, not frequency domain data. Filtering can also be implemented in the frequency domain (by multiplication), but the functions you're using don't do that.


How to compute impulse and frequency response of Flanger? I have to the implementation of a Flanger effect on Matlab, buy previously I have to plot its frequency response and impulse response.

The difference equation is $y[n]=x[n]+a\cdot x[n-d[n]]$
where $a$ is a constant, $|a|<1$, and $d[n]=\frac D2 (1-\cos(2\pi f_s n))$;
$D$ and $f \text{ const}$.

I'm having trouble on how to calculate any of DTFT o Z Transform of such difference equation. I can't find how to compute the varying time shifting transform for $x[n-d[n]]$.

Is there any other way to compute impulse and frequency responce?

Only linear time invariant systems are fully characterised by their impulse response. Modulation effects like a flanger break time invariance. While you can still create a momentary impulse response (which is only of very limited use), there is no meaningful frequency response.


Audible mp3 artifacts? What kind of compression artifacts would typically be most noticeable to a non-audiophile listener in a low-bit rate mp3 (MPEG-2 Layer III) compressed audio file?

What kind of audio signal might one synthesize that would make these artifacts maximally audible between "low" and "high" bit rate mp3 compressed files (using commonly available encoders)?

Pre-echo on very short and percussive sounds - for example crisp solo drums, or better, electronic drums. The TR-808 rimshot and clave are very short and impulsive sounds, about 5ms long (for the rimshot), and at less than 160kbps their attacks get "swallowed". At less than 96kbps some notes appear to be duplicated or out of place rhythmically because of pre-echo. Worse, each note in a repeated pattern will sound different depending on how it falls within the analysis window. Here is an example.

Try solo vocals too. I think casual listeners are better at detecting anomalies in voices (a kind of auditory "uncanny valley" effect). A classical example used in the audio coding community is the intro from Suzanne Vega's "Tom's Diner".


Variance of a histogram in matlab my doubt is about a vector of values that I have and which histogram is the following one:

The function I used to create the histogram was:

histogram(q,1000); %q is the aforementioned vector of values

Now, I have to calculate the variance of each of the 2 gaussian distributions that appear in the image. I know that the function  var(x) gives you the variance of a vector, so I tried with:

h = hist(q,1000);
h1 = h(1:500);
h2 = h(500:1000);

v1 = var(h1); v2 = var(h2);

But it's obvious that the variance of each one will depend of the vector I use to calculate it. Another possibility is that knowing the mean of each distribution (-1 and 1 respectively), I could evaluate the pdf of a gaussian distribution in $f_x(\mu)$ and I would obtain that:

$$f_x(\mu) = \dfrac{1}{\sigma \sqrt{2\pi}}$$

However, the distributions aren't normalised and the normalisation of each distribution would depend again of the vector I chose.

I have no more ideas, so I hope someone can help me.

Thanks in advance for responding.

Your distribution is bi-modal and not Gaussian. I'm assuming you want to find the variance of each of the two modes (which look Gaussian by themselves).

In your case the modes don't seem to overlap, so the easiest way would to split this into positive and negative values

q1 = q(q > 0);
q2 = q(q < 0);
v1 = var(q1);
v2 = var(q2);


Using Total Variation Denoising to Clean Accelerometer Data I know this is maybe a very basic question but I am doing this as a hobby and I can't find a solution to this problem.
Basically I am trying to remove some noise from data I am reading from an accelerometer. This is what I want to achieve (Taken from Total Variation Denoising (An MM algorithm)):

I read in Picking the correct filter for accelerometer data that Total Variaton Denoising would fit my needs. So I read Wikipedia - Total Variation Denoising article from Wikipedia and I think I have to use one of this equations:

But I don't understand how I apply this to my signal. Suppose I have a set of x,y points like in the plots above, how I apply the equation to that data? I implemented some simple low-pass and high-pass filters like this:

gravity[0] = alpha * gravity[0] + (1 - alpha) * event.values[0];

But this is maybe too complex and I don't know where to start or how. I want to implement this in Java or C so Matlab is not an option (I have seen a lot of MatLab implementing this). I will appreciate any help to guide me in the right direction!

Apart from Total Variation Denoising you could try a first much simpler approach: a median-filter. You just move a window along your data and replace the current input value by the median of all data in the window. You just have to optimize the window length (by experimenting).

By the way, the equations you copied into your question are for 2-dimensional data, but your data are 1-D.


bandwidth vs. bandwidth range As known, Bandwidth is the difference between the upper and lower frequencies in a continuous band of frequencies.
OK, what means when saying that bandwidth range is 15 ~ 20 KHz ?  Does that mean that the bandwidth should be 5KHz ?

Range is the minimum to the maximum value of something. Since bandwidth is the difference between the upper and lower frequencies in a continuous band of frequencies, then the range of bandwidth of 15 to 20 KHz would then mean that the minimum difference (bandwidth) would be 15 KHz and the maximum bandwidth would be 20 KHz. If the saying was instead that the range of occupied frequencies was 15 to 20 KHz then it would then mean the bandwidth was 5 KHz.


How to optimally normalize filter coefficients? I have designed an IIR filter in Matlab that comprises the following biquads
0.244   0.002   0.244   1   -1.432  0.543
0.244   -0.345  0.244   1   -1.555  0.777
0.244   -0.388  0.244   1   -1.652  0.943

I want to implement this filter in a 16-bit DSP. Before converting the coefficients to short int, I therefore normalize them by dividing all coefficients by the largest magnitude present in the matrix (in this case 1.652). This results in
b0      b1      b2      a0       a1     a2
0.147   0.001   0.147   0.605   -0.866  0.328
0.147   -0.208  0.147   0.605   -0.941  0.470
0.147   -0.235  0.147   0.605   -1      0.570

I then convert these coefficients to 16-bit values ranging from +32767 to -32768. I implemented the filter using the Transposed Direct Form II.

The I realized that the gain is very low although the filter shows 0 dB in the passband when testing it in Matlab.
The algorithm is as follows. It iterates through all biquads and reintroduces their ouputs to the next one as necessary.
for (k = 0; k < N_sos; k++) {

    y_n = (short) ((s1[k][1] + b[k][0] * x_n) >> 15);

    s1[k][0] = s2[k][1] + b[k][1] * x_n - a[k][1] * y_n;

    s2[k][0] = b[k][2] * x_n - a[k][2] * y_n;

    s1[k][1] = s1[k][0];

    s2[k][1] = s2[k][0];

    x_n = y_n;
}

IIR_out = y_n;

Notice that I have not taken the a0 coefficient into account. Since it usually is 1, it does not matter. But I wonder if the missing gain is caused by this missing coefficient. I tried multiplying y_n with it by exchanging
y_n = (short) ((s1[k][1] + b[k][0] * x_n) >> 15);

with
y_n = (short) a[k][0] * ((s1[k][1] + b[k][0] * x_n) >> 15);

but this heavily distorts the signal.

The OP should first implement the filter using floating point to confirm all specifications and structure. Once confirmed there is no "optimal" since the solution is a trade between precision required and accuracy needed in the result but typically the design approach is to use 2 more bits of precision for the coefficients over what is used in the datapath; scale all coefficients accordingly, let the filter grow the signal with an extended precision accumulator (such that no clipping occurs with input set to the maximum possible signal) and then scale the output back to the datapath precision.
See my further details on this at this answer recently posted:
What is the correct gain of an RRC Filter?


Correct method for processing FFT in realtime I am attempting to build a 'live' FFT for continuous analysis of an audio signal.  I am relatively new to this, but have attempted the full thing, and then updating the graph at 20Hz refresh rate.  My FFT is working with regards to identifying the correct frequencies, however is very 'jumpy' with regards to the amplitude. Before I analyse my code to see where I have gone wrong, I would like to just confirm the process I have used is correct.

Read audio samples into circular buffer.
Read half of the window size from last half of previously processed array.  (50% because I am trying to achieve 50% overlap).
Read half of the window size from the buffer to give the complete window size.
Perform Window Function on the window to process.
Zero pad the end of this to achieve a n^2 bit FFT Size.
Perform FFT on the full array
Save Array to List of arrays to average (add to end, delete first one to always retain x no. of averages).
Average the arrays in List of arrays (treat real and imaginary values as two independent values when averaging).
Obtain the magnitude for each bin in the averaged array by obtaining the square root of the imaginary^2 + real^2. This is then divided by the number of bins, and the entire value multiplied by 2 as we have halved the FFT spectrum to only look at the positive values. This is then multiplied by a window correction factor.
This thread should then sleep for however long it takes to capture half the window length needed (the next set of audio samples) subtract the computational time to perform this FFT step.

That is currently the method I am working to.  I know my FFT algorithm is good as I have tested it on known input and output values.  I am sure it is a coding mistake, but as this is my first project I wanted to check that the overall method I am following is sound.  If anyone can see any potential issues, please let me know.  Thanks.

There are numerous ways to do this and you have the general flow except, as noted, complex averaging is only done in special circumstance, which you haven’t mentioned and are in most cases not applicable.

One can average magnitude or magnitude squared. There is a 1975 Bell System Technical Journal article that goes into detail but a simple comparison is good enough for most people.

Zero padding a FFT to a power of 2 is probably unnecessary if you use a modern FFT library like fftw.  A library will also offer an optimized real FFT. The other advantage of not using a zero pad is that you can pick your bin widths more intelligently.  It may be convenient for a user to have something like 1Hz or 5Hz or 10Hz bin widths.

The number of averages should have some reasoning as well.

Depending on your actual application, some filtering prior to your FFTs might be useful. A seismologist might be interested in frequencies below 10Hz, but a lot of other people don’t. Some degree of pre whitening is often useful.

One can have a “correct” processing chain , but nevertheless underachieve as an analysis tool.

There is a lot about how to display this information that hasn’t been touched. Real time usually means real time display.

When I started doing DSP there were a lot of tricks that you can use to save cycles. You can avoid square roots by taking the max( real, imag) - \alpha min(real, imag) where you set any negatives to positives first.

Another trick is to complex bandshift filter to use a complex fft without redundancy


This is not an easy task since its a classification problem. A while back at work I used time-frequency-domain transformations like STFT (Short Time Fourier Transformation), IWT (Integral Wavelet Transformation), DWT (Discrete Wavelet Transformation) STCEP (Short Time Cepstrum), WVD (Wigner Ville Distribution) to create images/matrices from voice recordings which were then classified using a convolutional neural network as classifier. This is a common approach for this task. You could try that or you could try setting certain limits, e.g. if in a certain frequency range a certain energy is exceeded for a certain time, the file is classified as a speech file.


What is Causing a Cyclic 1kHz Noise Signal in Audio Recording? I dug my old Behringer UCA202 out of my shed and plugged it into a computer. I did some tests and it appears to be working fine... apart from some noise signals on the line input channels. I'm not sure where this type of noise could be coming from; weather it's the computer, the powered USB lead, or the device itself. What components could cause this type of interference?

###Spectrum Analysis and Statistics

###8 Second Spectrogram 750Hz - 5kHz (showing ~3.3Hz Period)

###1 Second 0-10Hz Low-Freq Noise Signal(related? i.e above ~3.3Hz cycle)

I have unplugged everything else from the computer (except power lead and VGA monitor).
I have tried plugging the device into a powered hub. The DC input on the hub states 5v 2.5A, but I'm using a 5v 2.0A power source. I don't know if this matters.
This is the recording:
https://soundcloud.com/n00dlz/uca202-testdevices-unplugged-monitor-off/s-H9zQsQzqArY
Only the monitor, mains lead and UCA202 are plugged into computer, and I unplugged the monitor from the mains while recording. There is nothing plugged into the UCA202.

"The UCA202 uses the ubiquitous TI PCM2902 integrated USB DAC chip"

So this a "Full speed" USB device, and the USB data is being sent at 12 MHz, but it doesn't send data continuously; it sends it in short bursts every 1 ms, which means it has a frame rate of 1 kHz.  This 1 kHz noise finds its way into everything and it's really annoying when you're designing USB audio products.

Sometimes it becomes much more audible because of a ground loop and you can improve it by lifting a ground somewhere, like unplugging your laptop from the wall, or connecting to a balanced input on your monitor instead of an unbalanced one, or using an audio transformer to isolate one thing from another.

  There is nothing plugged into the UCA202.

What if you short out the inputs instead of leaving them floating?  They're RCA inputs so they might not have normals that short themselves when not connected to anything.


How to correct a QAM time delay How can you correct the anticlockwise rotation that occurs in the received constellation of a QAM system?

My only idea is to compensate the delay with an advance but this is impossible since you would have to decode the message before receiving it.

If this represents the received signal where $a$ and $b$ are the real and imaginary parts respectively $s(t-\Delta) \approx [a \cos{(2\pi f_o \Delta)}-b \sin{(2 \pi f_o \Delta)}] \phi(t) \sqrt{2} \cos{(2 \pi f_o t)} + [b \cos{(2 \pi f_o \Delta)} + a \sin{(2 \pi f_o \Delta)}] \phi(t) \sqrt{2} \sin{(2 \pi f_o t)}$

one can see exactly how the signal has been rotated. Therefore the rotation can be countered by introducing a sine or cosine to counteract the delay.


How do you calculate the group delay in samples of a Biquad Filter at a given frequency from the coefficients? If I have a BiQuad all-pass filter with a given set of coefficients such as:

g: 0.903228
a1: -1.88142
a2: 0.903228
B1: -2.08299
B2: 1.10714

And my operation per sample to process the filter is:

    inputGained_2 = inputGained_1;
    inputGained_1 = inputGained;
    inputGained = input * g;
    output_2 = output_1;
    output_1 = output;

    output = inputGained + (B1 * inputGained_1) + (B2 * inputGained_2) - (a1 * output_1) - (a2 * output_2);

    return output;

And this is running at 44.1  kHz: What equation could I used to calculate the group delay at a given frequency like say 82 hz?

Thanks.

Pick two different frequencies slightly above and below your target frequency. The choice depends on what numerical precision you have available and what your sample rate is.
Calculate the z-transform at both frequencies
Calculate the phase at both points, make sure to check for "phase wrapping"
Subtract the phases and divide by the frequency difference
Scale properly and multiply with $-1$


How are CDMA codes distributed? I think this question is so trivial that it's nowhere explained (I would not find it) and the answer is probably trivial too. But I still wonder how codes in a typical CDMA system are distributed? They act like symmetric keys in a crypto system: both parties need to know them and there is a complexity of O(N) to share them (N being all possible users).

For example, in a WCDMA cellphone system, there is an "infinite" number N of possible nodes (e.g. a roaming user coming to the US). The codes can't be statically assigned. And the wireless base station can't just try all possible codes.

So I assume that the first communication of a handset and base station does not go over CDMA but the base station somehow detects that there is a new handset, creates a new code and sends this code to the handset. Afterwards the communication goes over CDMA.
For handover between different stations, the base station could distribute the code among its neighbor stations.

How does it work in practice?
And am I correct that the seed to an LSFR generating the chipping sequence is distributed as opposed to the chipping sequence itself?

I don't know how WCDMA does it, but CDMA2000 encodes the data twice or three times (it's been a while): once with a "long code" (42 or 43 bits, so it takes days to cycle through even at chip rates in the MHz), once with a "short code" (16 bits, with zero-stuffing so that the cycle is exactly 65536 chips long, and the spectrum is flat), and IIRC once with a 64-bit long Walsh code.  But There's a channel that's only encoded with the long code; that one is used by the base station to assign a specific set of codes and timings to the phone.

So, basically, it's done dynamically.

The CDMA2000 specs are out there.  When I looked (over ten years ago, which is why this answer is a bit fuzzy) the "real" spec was behind a paywall, but there were draft specifications that were pretty close.


Arbitrary Phase FIR design using FFT I want to design a digital FIR filter to cancel the phase of an analog Butterworth filter (i.e. the phase is smooth) i.e. the filter has unity gain with just an inverted phase component.

The following standard method for designing linear phase filters via FFT guarantees that the phase is linear (without deviation) even after applying a window (the window is obviously useful to prevent ripples occurring if the required impulse response is longer than the desired length):
FIR filter design with frequency sampling method - setting proper phase response for IDFT

Is there a similar method where by careful maths you can ensure that the magnitude is guaranteed to be unity (even after windowing to smooth any "ripples" in the resulting phase if the impulse is longer than the desired length)?

You are trying to design a FIR allpass filter. That's impossible (other than a simple delay) since true allpass filters consist of zero/pole pairs that are inverse of each others and the FIR poles are always at the origin.

However, with enough resources, you can always get "close enough" for your specific application. Of course, you would have to define what "close enough" means.

The actual design is simple enough

Sample your filter target on a very dense FFT grid
Do an inverse FFT
Window/truncate the resulting impulse response until you meet your "close enough" requirements. Typically the impulse response will decay exponentially, so truncation often gives better results that any type windowing (for the same filter length)
In most cases the impulse response will be non-causal, so you will need to add bulk delay.

Below is a piece of Matlab code that does that

%% start with a Butter filter
fs = 44100; % sample rate
bwCut = 1000; % 1000 Hz cut off frequency
bwOrder = 3;  % thrid order
%% design the filter
[b,a] = butter(bwOrder,bwCut/(fs/2));

%% sample on a large FFT grid
nFFT = 16384;
delta = zeros(nFFT,1); delta(1) = 1;
bwIR = filter(b,a,delta); % impulse response
bwTF = fft(bwIR);         % transer fundction

%% design the target
targetTF = exp(-1i.*angle(bwTF)); % that works every where except for Nyquist, where
                                  % the phase is undefined, we'll manually smootht this out
targetTF(nFFT/2+1) = 1;           % Nyquist phase repair
h0 = real(ifft(targetTF));        % target impulse response

%% now we have to cut it to a practical size
% let's see if we can capture all samples that are above a certain
% thershold. That's in essence a rectangular window
thresh = max(abs(h0)).*0.01;      % -40dB re max
tmp = circshift(h0,nFFT/2);       % shift to the middle of FFT window
i1 = find(abs(tmp) > thresh);     % find the area above threshold
h1 = tmp(min(i1):max(i1));        % grab that area

%% analyse the result. This really depends on the application.
% In this example we will look at the difference of the resulting
% transfer function to the "ideal" answer. The difference looks at both the
% real and imaginary parts, so it's a decent metric that combines magnitude
% and phase. It coould be interpreted as the amount of "noise" induced by
% the truncation error.
idealTF = abs(bwTF);

% model the results, rectangular window
h2 = conv(h1,bwIR);                   % convolution
h2 = h2(1:nFFT);                      % cut off the excess length
h2 = circshift(h2,min(i1)-nFFT/2-1);  % time align
diffTFRect = fft(h2)-idealTF;         % calculate difference
% same length hanning window
h1w = h1.*hanning(length(h1));
h2 = conv(h1w,bwIR);                  % convolution
h2 = h2(1:nFFT);                      % cut off the excess length
h2 = circshift(h2,min(i1)-nFFT/2-1);  % time align
diffTFHanning = fft(h2)-idealTF;

%% plot it
clf;
i = (1:nFFT/2);
fx = [diffTFRect diffTFHanning];
plot(i*fs/nFFT, 20*log10(abs(fx(1+i,:))));
xlabel('Frequency in Hz ->');
set(gca,'xlim',fs*[1/nFFT 1/2]);
set(gca,'xscale','log');
set(gca,'ylim',[-120, 10]);
ylabel('Difference in dB');
grid('on');
legend('Rectangle','Hanning');
title('Error');


Integer Wavelet with Only Positive coefficients Is there any Integer Wavelet Transform which produce only positive coefficients?
Currently i'm using Haar wavelet but it produces both positive and negative coefficients.

If you are talking about the wavelet by itself, no, because   it is a  zero-mean function (by construction), hence there are positive and negative values.

For wavelet coefficients from data, no, because of linearity: if for some $x$, $W(x)$ is positive, then $W(-x)=-W(x)$ is negative.

Note that for some signal models, with certain multiscale schemes, you could achieve your goal. You can find hereafter two references:

Lossless data hiding based on integer wavelet transform

  In data hiding, pieces of information represented by the data are
  hidden in the cover media. In some applications, people do care about
  the cover media. That is, the hidden data and the cover media may be
  closely related. For this type of data embedding, in addition to
  perceptual transparency, for some applications such as medical
  diagnosis and law enforcement, it is desired to invert the marked
  media back to the original cover media after the hidden data have been
  retrieved. The marking techniques satisfying this requirement are
  referred to as lossless, distortion-free or invertible data hiding
  techniques. From this point of view, it is observed that most of the
  current digital watermarking algorithms are not lossless.

Distortionless data hiding based on integer wavelet transform

  A novel distortionless image data hiding algorithm based on integer
  wavelet transform that can invert the stego-image into the original
  image without any distortion after the hidden data are extracted is
  proposed. This algorithm hides data into one (or more) middle
  bit-plane(s) of the integer wavelet transform coefficients in the
  middle and high frequency subbands. It can embed much more data
  compared with the existing distortionless data hiding techniques and
  satisfy the imperceptibility requirement. The image histogram
  modification is used to prevent greyscales from possible overflowing.
  Experimental results have demonstrated the validity of the algorithm.


How do I apply a function window to a signal? How can I apply a window function like Hamming or Lanczos to a signal, using its coefficients?

I mean, which method can I use to do this? FFT? Convolution? Which method has the better performance?

How can I apply a window function like Hamming or Lanczos to a signal, using its coefficients?

Just multiply, point-by-point.


How can I load a wav file, split it in frames, apply different filters do different frames and then recombine the filtered frames into a single file? I'm working on a wav audio file using MATLAB. What i would like to do is:

load a wav file and divide it into frames of the same dimension
filter these frames (with a different filter for each different frame) using, for instance, iirnotch
put all the filtered frames together again without having discontinuities between them

I already wrote a script using the overlap-and-add technique using a "Hann" window and it works quite well but between the different frames I'm not able to avoid the discontinuities so that when I listen to the wav file (made by adding these frames together again) I hear an annoying glitch disturbance.

Anyone know how to avoid this problem?

add technique using ann "hann" window and it works quite well but between the different frames I'm not able to avoid the discontinuities

That's probably because the filtering has tails that extend beyond the frame.  The window is supposed to taper the signal to 0 at the edges, but after filtering it is no longer 0 at the edges, which introduces discontinuities.

Are you zero-padding before filtering to keep these response tails?

I believe you need to apply a sqrt(window) before filtering ("analysis window"), and then another sqrt(window) after filtering ("synthesis window"), so the combination is a normal window, and the tails are tapered to 0, and then overlap the windows by the appropriate amount so that they sum back to unity.  Maybe also need to zero-pad before filtering so that the filtering tails don't wrap around, if you're doing FFT filtering, for instance.  More information here: http://www.dsprelated.com/dspbooks/sasp/Weighted_Overlap_Add.html

(If you're doing positive-phase time-domain filtering, it might be sufficient to use a regular window, zero-pad before filtering, and then instead of windowing again, just keep the entire zero-pad region with the tail in it when you sum back together, but I haven't tried it so I'm not sure.)


The Effect of a Kernel on an Image By Filtering I am going through an old exam paper and trying to understand this question where they ask us to tell them what effect the following image filter would have if its convolved with an image:

     0     0     0     0     1     0     0
     0     0     0     0     0     0     0
     0     0     0     0     0     0     0
     0     0     0     0     0     0     0
     0     0     0     0     0     0     1
     0     0     0     0     0     0     0
     1     0     0     0     0     0     0

Now by using Matlab I can easily tell that this blurs the image. But what I am wondering is how would one know by just looking at this filter that the image will get blurred when this filter is applied?

Here's the image I'm working with in its original form:

Here's the same image after being convolved with the kernel in Matlab:

I would appreciate an explanation here or if someone can point me in the right direction as to a reading online I'd appreciate that too.

Convolution with following kernel returns the signal unchanged.

 0     0     0     0     0     0     0
 0     0     0     0     0     0     0
 0     0     0     0     0     0     0
 0     0     0     1     0     0     0
 0     0     0     0     0     0     0
 0     0     0     0     0     0     0
 0     0     0     0     0     0     0

If the '1' is shifted to the right one position, the convolution output will be shifted to the right one position as well.

If the '1' is a '2' the convolution output will be double (scaled by 2).

If there is more than one non-zero number in the kernel, the convolution output will be the addition of each scaled and shifted copy.

For your kernel

 0     0     0     0     1     0     0
 0     0     0     0     0     0     0
 0     0     0     0     0     0     0
 0     0     0     0     0     0     0
 0     0     0     0     0     0     1
 0     0     0     0     0     0     0
 1     0     0     0     0     0     0

the output is the addition of three copies of the signal, shifted by (1,3), (3,-1) and (-3,-3) and scaled by 1. I wouldn't call it blurring though.


method to detect frequency sweep in audio I have search a bit to try to find a similar questions, but I have not been able to find any leads.

There are techniques to programmatically detect steady tones in audio, such as the Goertzel algorithm. I have been trying to find a way to detect a frequency sweep of some kind. My specific application is to monitor the audio from a security camera, and detect, say, when a door prop or car alarm is triggered. If it matters, one can assume that the sound will repeat many times, and that the frequency is a smooth sweep, and not discreet tones.

I don't want to base it strictly off of amplitude, since the frequency response is sometimes very non-linear on these cheap cameras, and other noises can be a source of false alarms. Also, the alarm source may be distant.

If you know the signal of interest in advance, you can pre-record it and constantly calculate correlation between recorded sample and incoming signal.

Record a sample sound from a channel to be analyzed (you can record a number of periods)
Reverse it
Use reversed signal as an impulse response for FIR filter (this also is known as "matched filter") -  filter will output a peak when signal of interest passes through it.

If you are interested in general varying tone detection (much more difficult problem), look at the chirplet transform - the sort of wavelet transform which uses varying frequency wavelets (chirplets) as its basis. People used it to count wolf howls, for example.


What is the best Algorithm to estimate the frequency for Almost-Periodic Data or complex periodic data? I’ve got what it looks like a periodic signal but the periods seem not equally repeated, I want to estimate the PSD, particularly I need to obtain the dominant frequency of this data. I’ve estimated the autocorrelation and then applied FFT, the results as shown in the figures below.
It seems that FFT algorithm is not the best option for this kind of data. Could anyone please suggest another solution?.

Time series

autocorrelation

PSD

The dominant frequency always have the highest power and by looking at the PSD plot you can easily compute it. The dominant frequency in your case is where the PSD has the highest peak. In this case it is somewhere at around f = 2.7.


Where do we use windowing? I have been studying about the box window and am unable to understand one thing, do we first sample a signal and then window it resulting in up-scaling or downscaling, or is windowing a practical case of sampling itself which is performed on a continuous time signal?

Alright, the way to think about the rectangular window (as well as all of the other windows) is that you begin with an infinite sequence of samples and you multiply the samples outside the window with 0 and the samples within the window with something non-zero (and usually, but not always, positive).

So the DFT (of which the FFT is fast method of doing the DFT) does periodically extend the $N$ samples passed to it.  But there are a couple of different ways of looking at it.  One way to look at it is with just the definitions of the DFT and iDFT:

$$ X[k] = \sum\limits_{n=0}^{N-1} x[n] \ e^{-j 2 \pi nk/N}  $$

$$ x[n] = \frac{1}{N}\sum\limits_{k=0}^{N-1} X[k] \ e^{j 2 \pi nk/N}  $$

The $x[n]$ that the DFT $X[k]$ is aware of is simply a periodic sequence with period of $N$. both $$ x[n+N]=x[n] $$ and $$ X[k+N]=X[k] $$ are true.

Now the alternative way to think about the DFT is that it is the sampling of the spectrum of the Discrete-Time Fourier Transform (DTFT):

$$ \hat{X}\left( e^{j \omega} \right) = \sum\limits_{n=-\infty}^{+\infty} \hat{x}[n] \ e^{-j \omega n}  $$

where $ \hat{x}[n] =
\begin{cases}
0, & n < 0 \\
x[n], & \text{for }0 \le n < N \\
0, & N \le n
\end{cases} $

and this is the same as multiplying by a rectangular window.

$$ \hat{x}[n] \ = \ x[n] \ w[n]  $$

where $ w[n] =
\begin{cases}
0, & n < 0 \\
1, & \text{for }0 \le n < N \quad \quad \\
0, & N \le n
\end{cases} $  (length-$N$ rectangular window)

But here you cannot tell that you extracted $\hat{x}[n]$ from a periodic $x[n]$ or if those $N$ samples came from an infinitely long and non-repeating stream of samples.  Whatever it is, $\hat{x}[n]$ is non-repeating, infinite in length (with all those zeros padded), and you can apply the DTFT to it.

It doesn't have to be the rectangular window, it could be any length-$N$, finite-valued window for $w[n]$ and the math that follows is the same.

Now with that definition of $\hat{x}[n]$ and  $\hat{X}\left( e^{j \omega} \right)$, it's not hard to show that

$$ \hat{X}\left( e^{j \omega} \right) \Bigg|_{\omega = 2 \pi \frac{k}{N}} = \hat{X}\left( e^{j 2 \pi k/N} \right) = X[k] $$

Now, what happens in either domain, is that when you uniformly sample the function (say, in this case, the spectrum $\hat{X}\left( e^{j \omega} \right)$)  in one domain, it causes in the other domain, the signal ($\hat{x}[n]$) to be repeatedly shifted (by $N$ samples) and overlap-added:

$$ x[n] = \sum\limits_{m=-\infty}^{+\infty} \hat{x}[n+mN] $$

It is here where the discontinuity between $x[N-1]$ and $x[N]$ (which is the same as $x[0]$) becomes a concern, and tapered windowing to reduce or eliminate this discontinuity is perhaps indicated.


Slow Down Music Playing While Maintaining Frequency Playing a piece of music audio at a slower speed would lower its pitch (frequency). Is there a tool and theory to slow down the song playing while keep the frequency the same? I suppose one can do windowed Fourier transform or wavelet transform. It seems one have to either pre-select the window size or dynamically choose the wavelet basis. Is there any specific and detailed theory and application for doing that?

Yeah some of us can do it, you can speed up or slow down without affect the pitch, some guys call this applications of Time Stretch, there different ways to do it, you can do in frequency domain or time domain, you will need choose what is best for you, you will find some advantages and disadvantages of each.

Time Domain:

In Time Domain you can try some techniques like:

TDHS(Time domain Harmonic Sacaling)
SOLA (Synchronous Overlap Add)
PSOLA(Pitch Synchronous Overlap Add)
WSOLA(Waveform Similarity Overlap Add)

Pros: Is fast, some algorithms are easy to understand, good quality in monophonic sounds.

Cons: Generally you'll need a very nice pitch track to splice in the right position, it is hard to do :-(, so if your pitch track fails or not work in Poliphonic sounds this algorithms will give a lot of glitches/artifacts in the output sound.

Frequency Domain:

All time stretch that i know in frequency domain are based in phase vocoder techniques.

Pros: Will work in polyphonic or monophonic sounds.

Cons: Can be painful understand all the math, implementation is a bit hard, is not so fast like time domain codes, for voice i prefer time domain results, some tricks to improve the result of the standard phase vocoder are not shared.

I can say that the window and the hop size are one of the key for the phase vocoder quality, generally we choose 4x overlap to resynthesis, one hann window of size 4096 is enough for my ears (of course if u have processing power for this sizes), the standard phase vocoder can add some reverberant effetcts, to try avoid this kind of problems you maybe need lock the phase.

For datailed take a look in the Miller Puckette and Portnoff paper


Random (Over) Sampling signal and perfect Reconstruction in Nyquist form? Imagine we have band limited signal with bandwidth of $B$, so the required Nyquist rate would be $f_{nyq}>2B$ that is oversampled with rates $f_s$ where $f_s = M*f_{nyq}$ and $M$ is random and $M>=1$. By that, I mean the sampling rate is at least at Nyquist rate but most of the time its way more than Nyquist and it is random. I have heard and read about  iterative  My questions are:

Is there a method to construct the signal as if it was taken in Nyquist more, without degradation of signal quality? After all, we have more samples than needed, so I do not expect any degradation, right?
Is there are some loss, how much is it, what factors it depends on?
Is there a toolbox or something that I can feed my data in and have some experiments of my own? ,

4.Is the technique used in MATLAB lossy?

I appreciate any kind of contribution here :)

For the moment restrict the discussion to 1D signals. Sampling operation can be performed as uniformly or non-uniformly. The uniform sampling is the most obvious, simplest and the preferred form, unless otherwise stated.

Nonuniform sampling can be performed in a number of ways, periodically or aperiodically, sampling instant given by a formula or being completely random.

Depending on the nonuniform sampling form, the reconstruction criterions may change. For example when you have a formula based nonuniform sampling instant $t_n$, you can (if possible) find a kernel $\phi_n(t)$ for the reconstruction integral that provides exact (perfect) recovery of a continuous time signal $x(t)$ from its nonuniform samples $x_n$. In this context the $sinc()$ kernel is a special case of the generalized sampling theorem WKS (Whittaker-Kotelnikov-Shannon) with sampling instants given by $t = n T_s$. This was historically called as the Nyquist sampling, who most successfully emplyed it into the earliest digital pulse communication systems.

The primary requirement that Nyquist (uniform) sampling theorem would provide exact reconstruction of a bandlimited baseband signal $x(t)$ from its uniform samples is that the sampling rate satisfy the Nyquist criterion; $F_s \geq 2 F_c$ where $F_c$ is the signal's bandwidth.

Now for the non-uniform sampling strategies, the most common restriction for perfect signal recovery is that on the average the Nyquist rate associated with the bandlimited signal is maintained. So when there is a deficiency of samples at some short interval $I_1$ then there must be an excess of samples at another closeby short interval $I_2$ so that the average Nyquist rate preserved.

The above condition can be stated as the local Nyquist rate criterion.

Note that unlike uniform sampling case, for nonuniform sampling strategies and exlcusively true for the random sampling strategy, the sampling times $t_n$ in addition to the samples $x_n$ should be known in order to reconstruct the signal perfectly.

You can look for Whittaker's, Yen's, Shannon's papers for classical view on the nonuniform sampling of signals.


What does SDMA's space actually mean? What does SDMA's space mean?
I mean,T for TDMA is time,and unity of time is second. F for FDMA is frequency,and unity of frequency is Hz.Now i know the S for SDMA is space or spatial,but i don't actually understand what does this space mean?is it mean the area,so S's unity is $m^2$ ?

It means that in a multipath fading environment, the signal will have several dominant spatial paths from transmitter to receiver.

From the perspective of a receiver, this means that the received signal will have components from several angles of arrival (AoA).

MIMO techniques exploit this, and use joint beamforming between arrays at both the transmitter and receiver to send different data streams along these spatial paths.  To get the full multiplexing gain, however, these signal paths must be statistically independent.

If the environment does not have any multipath characteristics (i.e. free space), then it is not possible to provide spatial multiplexing, and the best you can do with multiple antennas is conventional beamforming (CBF).


Guitar effects explanation books Can you please suggest some books to get how guitar effects: distortion, overdrive, phaser and so on are achieved ?

DAFX is a pretty good one.

Harmony Central used to have a surprisingly good knowledgeable base for algorithms but that seem to have disappeared since Guitar Center took them over.

A lot of source code can be found at http://www.musicdsp.org/, but it's fairly varied in quality.


How can I plot the spectrum of a signal in MATLAB? Consider an input signal x_total and the complex envelope of an output signal y_total in MATLAB:

T = 10e-9;
f_s = 1e12;
t = [0 : 1/f_s : T]';
f_c = 11.5e9;

f_1 = -0.1e9;
f_2 = 0;
f   = [f_1 f_2]';

x_1 = 1*sin(2*pi*(f_c + f(1))*t);
x_2 = 1*sin(2*pi*(f_c + f(2))*t);
x_total = x_1 + x_2;

y_envelope = A_out .* exp(1i* theta) .* exp(1i*Phi);

2 Questions:

How do I plot the spectrum of x_total? X_total = fft(x_total) yields a complex result. I want the x-axis of my spectrum to be in Hertz.

And I want to plot the spectrum of my output signal y_total as well. From the complex envelope, can I say y_total = A_out*sin(2*pi*f_c*t + (theta+Phi))?

I just read:
A = abs(X_total) * 2/(f_s*T) should yield the amplitude
and f = linspace(0,f_s,length(X_total)) should yield the corresponding
frequency.
Why do I have to normalize the amplitude by 2/(f_s*T) ? Applying this code yields the spectrum for x_total, which I want.
Applying the same for y_total does not yield a spectrum at all, when I plot it. Why not?

One way is simply to plot the absolute value:

X_total = fft(x_total);
figure;
plot(10*log10(abs(X_total).^2));

Note that if your signal is real, you can plot only the first part of the vector.
Another option is to use the Welch Periodogram

figure;
pwelch(x_total);


Causality as applied to capacitors This question stems from a point of confusion that I still have about the causality, linearity, and time-invariance in LCCDEs. I wanted to use the capacitor as an example.

Consider a capacitor with capacitance $C$. Taking the current $i(t)$ to be the input to the system and the voltage $v(t)$ to be the output we have $$i(t) = C \frac{\mathrm{d} v(t)}{\mathrm{d}t}$$

This differential equation can be solved to obtain
 $$v(t) = v(t_0) + \frac{1}{C} \int_{t_0}^{t} i(\tau) \mathrm{d} \tau$$

My first question is: isn't this mathematically valid for all $t$? In other words, does this give us the response for all $t$ or is it only valid for $t > t_0$? If it is valid for all $t$, including the $t < t_0$ case, doesn't this make the system non-causal since it anticipates future input and output values? Are we allowed to integrate backwards in time?

My second question relates to the assertion that a for the LCCDE to describe a linear system, the initial conditions must be zero. Suppose $t_0 = 0$ such that  $$v(t) = v(0) + \frac{1}{C} \int_{0}^{t} i(\tau) \mathrm{d} \tau$$

With $v(0) = 0$ the system is linear. But the choice of $t_0 = 0$ is arbitrary, since for example
$$v(t) = v(0) + \frac{1}{C} \int_{0}^{t} i(\tau) \mathrm{d} \tau = v(2) + \frac{1}{C} \int_{2}^{t} i(\tau) \mathrm{d} \tau$$

Why shouldn't we require that $v(2) = 0$ as well for that matter? What am I missing here? Thank you in advance.

Your first equation is indeed valid for all $t$, but this has nothing to do with causality. You have

$$\begin{align}v(t)&=\frac{1}{C}\int_{-\infty}^ti(\tau)d\tau\\&=\frac{1}{C}\int_{-\infty}^{t_0}i(\tau)d\tau+\frac{1}{C}\int_{t_0}^ti(\tau)d\tau\\&=\frac{1}{C}\int_{-\infty}^{t_0}i(\tau)d\tau-\frac{1}{C}\int_{t}^{t_0}i(\tau)d\tau\end{align}$$

If $t_0>t$ you simply subtract from the first integral the contribution from the interval $[t,t_0]$. But this is still causal, you don't need the future input, you just need to subtract it in case you added it in the first integral. The output could of course be computed without it.

Concerning your second question, you choose the time $t_0$ as the beginning of your processing. The value of $t_0$ doesn't matter but once you choose it, it's fixed. The system is linear if the initial conditions at $t=t_0$ are zero.


Techniques Used To Process Accelerometer Recordings I'm currently working on a project which involves using Accelerometer data. Having searched the internet for techniques in processing the data, I am unable to find simple techniques. Is there any suggestions, of what techniques to use in processing recordings? The data recorded is to do with body motion, walking up-stairs, down-stairs, falling etc.

The advantages and disadvantages of techniques to processing the data would also be welcome.

One suggestion (non-simple):  Obtain thousands of recordings from a broad spectrum of thousands of people engaged in your chosen activities to measure (with the activity externally logged, witnessed log book, video, etc.).  Feed all that data to several contemporary machine learning algorithms (DNN, etc.), and hope for a resulting inference method that meets your targets.


Given a discrete time signal, what is the sequence of possible frequencies I can get from DTFT? I know that when I have a discrete time signal, let's say:

$$f(t_n),\;t_n=\dfrac{n}{F_s}$$

The definition of the DTFT is given by:

$$F_n(\omega)=\sum_{n=0}^{N-1}f(t_n)\cdot e^{-i\omega_nt_n}$$

Now, my question is regarding $\omega_n$.

I know the frequencies will be discret because we can't measure lower frequencies than the difference between $t(n+1)$ and $t(n)$.

However, I am not clear how exactly this works, and how can one determine the "frequency vector". I am not sure I got the expression for $t(n)$ correct.

I am sorry this question isn't 100% clear, it's because I'm trying to make heads and tails from this definition and articles online are not very clear.

Thank you.

The discrete-time Fourier transform (DTFT) is defined by

$$F(\omega)=\sum_{n=-\infty}^{\infty}f[n]e^{-jn\omega}\tag{1}$$

Note that this is the common definition, and it is different from the one in your question. First of all, the sequence $f[n]$ doesn't need to be of finite length. Second, the frequency variable $\omega$ is a continuous variable. Your formula can be derived from $(1)$ by assuming that $f[n]$ has a finite length, and by computing the DTFT on a finite set of discrete (angular) frequencies. (Note that the frequency variable $\omega$ in $(1)$ is normalized by the sampling frequency, unlike the discrete frequencies $\omega_n$ in your question.)

Another related transform is the discrete Fourier transform (DFT):

$$\tilde{F}[k]=\sum_{n=0}^{N-1}f[n]e^{-j2\pi kn/N},\qquad 0\le k<N\tag{2}$$

which, for finite length sequences of length $N$, is just a sampled version of the DTFT:

$$\tilde{F}[k]=F\left(\omega_k\right),\qquad \omega_k=\frac{2\pi k}{N}$$

Unlike the DTFT, the DFT can be computed efficiently. Efficient algorithms for computing the DFT are referred to as Fast Fourier transform (FFT) algorithms.

If this doesn't clear up your confusion, leave a comment below.


Does instability make an otherwise LTI system nonlinear (or time-variant)? I am spinning this question off from the question from johnny.  Matt L. and I have had directly opposite conclusions to johnny's question.

I want to decouple the question from issues of causality and other goofy stuff.

So we have a simple first-order recursive system described with time-domain I/O equation:

$$ y[n] = p \cdot y[n-1] \ + \ x[n] \quad \quad \forall n \in \mathbb{Z} $$

Of course, the Z-transform of this is

$$ Y(z) = p \cdot z^{-1} Y(z) \ + \ X(z) $$

and transfer function

$$ H(z) \triangleq \frac{Y(z)}{X(z)} = \frac{z}{z-p} $$

We would normally identify this as a simple and realizable LTI system with a zero at $0$ and a pole at $p$.  But in the other question, there is an issue regarding linearity and time-invariance for the case when $p=-1 \ $.

For what values $p$ is this system linear?  For what values $p$ is this system time-invariant?

This is, I believe, the kernel of the disagreement I have with Dr. Matt L.

I don't know how the discussion arrived to this point but it is pretty convoluted to follow after this much of time. It is an ordinary difference equation with constant coefficients thus defines a linear, time invariant system. There is no need to pursue further from that if a unique solution exists. Here we have a problem with the emphasized part.

Let's first write the descriptor state space equations: the system is described by

$$
\begin{align}
\begin{bmatrix}1&-1\\0&0\end{bmatrix}
\begin{bmatrix}s_1[k+1]\\s_2[k+1]\end{bmatrix} &=
\begin{bmatrix}p&0\\0&0\end{bmatrix}
\begin{bmatrix}s_1[k]\\s_2[k]\end{bmatrix} +
\begin{bmatrix}0\\1\end{bmatrix}x[k]\\
y[k] &= \begin{bmatrix}1&0\end{bmatrix}
\begin{bmatrix}s_1[k]\\s_2[k]\end{bmatrix}
+[1]x[k]
\end{align}
$$

Now this is where I think the problematic part is. This system, though being LTI, is not regular (relevant buzz words are regular, impulse-free, index-1 descriptor systems). In other words, there exists no $\lambda$ for which the expression $\det(\lambda E-A)$ is nonzero and hence one of the modes is $\frac{\bullet}{0}$ and in fact this is $\frac{0}{0}$ for this example. This means that our system has solution uniqueness issues, unlike the causal LTI system, there is no guarantee on the existence of a unique solution. There is no guarantee of an admissible solution for that matter(buzzword impulsive systems). Hence the LTI reasoning of the other answers won't cut it.

And this what causes the trouble as far as I can tell from Matt L's argument is that he found two nontrivially distinct solutions for the same system and concluded that this cannot be a linear system. But this also assumes the uniqueness and existence of a solution and initial conditions.

It only differs from the regular systems in the way that the uniqueness and existence guarantees of the standard LTI systems cannot be assumed. The models no longer can be assumed to have admissible trajectories for all possible signals.


T wave analyis and extraction using discrete wavelets I'm doing a project in T-wave extraction from ECG signal. Initially I completed the preprocessing steps like noise removal and baseline drift cancellation. Then I need to extract T wave using discrete wavelet transform. I tried different wavelets for different level in wavelet, but I'm not able to extract the T wave.

Can you help me with MATLAB coding?

First, this is not an easy task to automate. An algorithm that works great on clean ECGs can fail miserably on noisy ones. That said, I would suggest you look at the approach used to identify the ST segments in the Long Term ST Database. It is one of, if not the, largest ST segment database publicly available. They describe the approach they used, which includes automatic segmentation followed by expert review.

Also, take a look at the PhysioToolkit from PhysioNet. I'm pretty sure the have some automatic ECG annotation tools. I looked, but couldn't find them. Still, I think they are there somewhere.


Why does sampling frequency affect the filter results? I have a big issue,  basically my goal is to create a low-pass butterworth filter  around 0-35Hz  with a sampling frequency of 5000Hz.
My Matlab code is

 Fs=5000;
Wp=30/(Fs/2);
Ws=45/(Fs/2);
Rp_db=-20*log10(.95);
Rs_db=-20*log10(.05);
[order,wn] = buttord(Wp,Ws,Rp_db,Rs_db);

But it doesn't work, becouse the filter explodes

As you can see the filter is unstable.But if i change the sampling frequency from 5000Hz to 500Hz things change

with the cutt-off frequency 0.135 (33,75 Hz).

Why does the sampling frequency affect so much the filter result?
Thanks in advance

Limited numerical precision. The higher the sample rate, the closer the poles move to the unit circle, the closer to the unit circle, the less stable the filter is.

There are different implementation methods that are better than others: design as poles, and zeros and not as transfer function, use cascaded second order sections, use correct section ordering, higher precision data types, etc.

In most cases very high order filters with poles close to the unit circle are a bad idea. You get very high group delay and phase distortions and massive time domain ringing.


Mat class of c++ in open cv I am a beginner in opencv and have started learning it using visual studio 2012 and opencv 2.4.7 . But i got stuck at very beginning .I want to know about Mat used in Opencv .It is a class and an image name is a variable of this class .i just know this .but still i am feeling like i didn't get it totally like what are its data members and member functions?

Mat declares a variable of type Matrix which can have multiple dimensions. In VS

Mat image = imread('test.bmp'); //lets consider a 2d image for now

Now put a . in front of image and the member functions will be exposed.


Reconstruction of BandPass Signal from Complex Envelope From what little I've been able to gather, the complex envelope preserves all information required for the reconstruction of the original signal and that it basically sends a low pass version of a bandpass signal.(Correct me if I'm wrong).I need an explanation for how the complex envelope is being created as per this text book? And also of how the original signal is obtained at the end?

Let’s get one thing straight. You cannot recover the original real and imaginary parts of a complex signal based only on the magnitude values of the complex signal. That’s equivalent to saying, “I have the magnitude of a complex number. How do I compute that complex number’s real and imaginary parts?” There are an infinite number of possible correct answers to that question!

However, as MBaz rightly said, you don’t have any “envelope” signal in your block diagrams. And by the way, there are no such things as “complex envelope” signals. All “envelope” signals are real-only. Having said all of that, you can learn everything about your two block diagrams by writing down the algebraic equations for all their various signals. (Your textbook should have done that for you.) Then pay attention to all the various trigonometric identities in your math reference book. It’s a mild pain to do all of that, but it’s terrifically educational.

And another thing, your two block diagrams are not the normal “complex modulator” and “complex demodulator” used in modern digital communications systems. In a normal digital comms system the real-valued gI(t)/2 and gQ(t)/2 would be added (or subtracted, I forget) to produce a new real-valued signal that’s transmitted via an antenna. Have a look at Figure 1 of the following PDF file: http://www.testequity.com/documents/pdf/keysight/complex-modulation-generation-wp.pdf


What do troughs represent in a sound wave? What's the notion of negative amplitude? I hope the title and the question don't mean different things.

Sound is a pressure wave in a fluid medium. That medium has a "static"or average pressure and sound is simply a small change of that pressure (as a function of time). Negative sound pressure just means that the current pressure is a little lower than the average pressure.
For example, static air pressure is around 100,000 Pascals (Pa) . Normal audible sound is quite a bit smaller, maybe 0.1 Pa (for 76 dB SPL). That means that the pressure are the microphone (or your ear) varies from 99,999.1 Pa to 100,000.1 Pa


Which digital modulation scheme would be suitable for transmitting digital data through the output of a PC soundcard I need to output a data stream at $40\textrm{ kbit/s}$ or faster out of a computer system whose only accessible output peripheral is the audio interface. This interface has reasonable specs, $96\textrm{ kHz}$ sample rate, with $24$-$\textrm{bit}$ resolution, but the output stage is AC-coupled. Only one output channel is available. A good assumption would be that it behaves like a band-pass filter with a $4\textrm{ Hz}$ to $40\textrm{ kHz}$ pass-band with less than $-1\textrm{ dB}$ attenuation; and that it has a $90\textrm{ dB}$ SNR. I don't have any other complexity constraint on the emitter.

I do not expect additional noise / attenuation on the cable connecting the emitter to the receiver.

The receiver is an embedded system with a $120\textrm{ MHz}$ Cortex-M3 MCU. If necessary, similar audio acquisition performance can be assumed. An additional dedicated demodulation chip (if there exists such a thing for such low frequencies) could be an option.

Which digital modulation scheme would be suitable for this situation?
Are there already code libraries (software defined radio libs?) that would prevent me from reinventing the wheel?
Are there existing applications with similar constraints at which I could look for inspiration?

You've got a pretty good set of circumstances here; you should be able to meet your goal without too much trouble. I don't see anything in your description that would eliminate a whole class of modulation (e.g. phase-shift keying, frequency-shift keying, etc.). Some of the factors that would go into the choice of a suitable format would include:

The required spectral efficiency (i.e. how much data throughput you need relative to the available bandwidth)
The complexity requirements for your receiver (which is usually the most complicated part of the system)
How much effort you're willing to put into developing the implementation.
Other circumstances specific to your application (e.g. if you have poor timing accuracy at one or both ends, known interference, or a poor channel response)

So, ticking these off one by one for your system, we can come up with a few guidelines:

It seems that your biggest constraint is your channel response (which is limited by your sound card's DAC). If you have 40 kHz of one-sided bandwidth available, then you're going to be limited to a symbol rate that is somewhat below that. For a target data rate of at least 40 kilobits per second, you'll want some scheme that transmits multiple bits per symbol.

Provided that your embedded platform isn't loaded down with too many other functions, a modern 120-MHz ARM processor should easily be able to handle demodulation of most any format in the tens-of-kilobits-per-second range.
I'm not sure specifically which model you're running with, but many recent processors provide very tight integration of onboard ADCs with the memory and interrupt subsystem, perhaps allowing you to (without manual CPU intervention) automatically sample the input signal at a specified rate, store the samples in onboard memory, and trigger a processor interrupt only when a certain-sized block of samples is available to be processed. I know that some Atmel devices at a minimum provide this sort of functionality; I've had good success with them in the past.

I assume that this is probably one component of some larger system, so you don't want to have to build something really complex in order to meet your requirements. Luckily, based upon what you've said, I think you should be able to implement something that meets your requirements pretty simply. I assume that since the two devices are cabled together, you're likely to have a very high signal-to-noise ratio (40 dB or greater), which for digital modulations is typically measured as $\frac{E_b}{N_0}$. That makes life much easier, as you likely could skip additions like error-correction encoding altogether (or at least implement a very simple scheme that would catch the very infrequent errors that you would expect at such an SNR).

As far as special circumstances go, for this system, I wouldn't expect much of any. I would expect the oscillator accuracy on the PC side to be pretty good (at a minimum crystal-controlled, so you're in the range of <50 ppm or so; possibly much better if the oscillator is calibrated using some other more precise source). The embedded side is likely to be the same; I assume you're using a crystal oscillator as the clock source. Since the two ends are cabled together, I'll assume that you don't have interference of note.

So, smashing all of these together into a single recommendation, I would probably start down the path of a quadrature phase-shift keying (QPSK) approach at 24 kilosymbols per second. At 2 bits per symbol, this yields a data rate of 48 kilobits per second, which is in excess of your requirement. This particular rate makes your implementation a bit easier; since the output DAC is running at 96 kHz, this results in 4 samples per symbol (it's always easier to run at an integer number of samples per symbol time). I would probably try to design the embedded side so that it would sample at the same 96 kHz rate if possible; this avoids the need to do any resampling on the more resource-starved end.
To avoid any problems with the DC notch that your sound card DAC employs, you could modulate the QPSK signal onto a carrier at 24 kHz. Then, the modulated signal's spectrum would have a null at DC, which would line up with your notch. It's possible that the notch could end up not being a problem at all (especially if it is truly just a few Hz wide as your suggested). In that case, you could potentially get by with an even simpler scheme that just works at baseband, bypassing carrier modulation altogether.
QPSK is a good choice because of its simplicity, both at the transmitter and receiver. At your SNR, you could accomplish greater spectral efficiency using a more complicated scheme like quadrature amplitude modulation (QAM), but the constant-envelope property of PSK signals is attractive from a receiver complexity standpoint. As a note, if you truly needed more bits per symbol in the future, you could move to a higher-order PSK constellation like 8- or 16-PSK. However, these are suboptimal from a bit-error-rate performance standpoint when compared to QAM constellations.
As far as a library implementation goes, I'm not aware of anything that you would just be able to drop in and go, especially for an embedded platform. Your receiver implementation is likely to be tied in with the hardware interface to some extent. You may be able to find some existing implementations for the various steps required for the demodulator, but you'll need to at least tweak what you might find to work well on your platform. The GNU Radio project is a good place to look if you just want to see C++ implementations of many different communications signal processing operations, and it might even yield a useful framework for implementing the transmitter onboard your PC. As a summary, the high-level steps that your receiver would need to do would include:

If a nonzero carrier frequency was used:

Frequency synchronization: locate and track the carrier frequency offset due to oscillator mismatches between the transmitter and receiver (in many cases, that frequency offset will be approximately constant over time)
Carrier demodulation: translate the signal to baseband, yielding two in-phase and quadrature (I/Q) signals (often expressed as a complex baseband signal).

Matched filtering: pass the baseband signal through a filter matched to the pulse shape used at the transmitter (you can likely get away with a rectangular pulse)

Timing synchronization: locate and track the times that correspond to symbol transitions; this can be done by tracking the peak locations in the matched filter output modulo the symbol time

Bit slicing: convert the matched filter outputs at symbol sample times to hard bit decisions

Serialization: make sure you write the multiple bits per symbol out in the right order!

This might sound like a complicated process, but building a practical receiver for even a simple situation like this can be very illuminating. Just comment if there's anything else I left out.


The spectrum of the sum of two periodic signals Considering the continue, periodic signals $s_1(t)$ and $s_2(t)$, with the period $P_1$ and $P_2$ respectively.

Consider now a new signal $s$ the sum of the two pervious signals.

I come to know that the spectrum of the sum, is not necessarily the sum of the spectrums unless the two signals $s_1(t)$ and $s_2(t)$ have the same period. With spectrum I mean:

$$X\in [\mathbb{Z} \rightarrow \mathbb{C} ]$$

$$X[k] = \langle x(t), \phi_k\rangle = \frac{1}{P} \int_P x(t)e^{-ik\omega_0 t}$$

with $\omega_0 = \frac{2\pi}{P}$. $P$ here represents the period.

Note that $x(t)\in [\mathbb{R} \rightarrow \mathbb{C}]$ is periodic with the period $P$.

I don't understand why? Because at it seems this violates the linearity properties.

the Fourier Series for a single periodic function, having period $P$,

$$ x(t) = x(t+P) \qquad \forall \ -\infty < t < +\infty $$

is

$$ x(t) = \sum\limits_{k=-\infty}^{\infty} c_k \, e^{j(2\pi k/P)t} $$

each harmonic $c_k$ is labelled with its harmonic number, $k$, in the above. the actual frequency of each harmonic, in Hz or units of 1/time, is $\frac{k}{P}$.  the fundamental frequency is the reciprocal $\frac{1}{P}$ of the period.  in angular frequency (radians per unit time), multiply both by $2\pi$.  but, in the harmonic number $k$ there lacks enough specific information to tell you what the actual frequency is of the harmonic.  e.g. the second harmonic of a periodic tone of 220 Hz is the same frequency as the fourth harmonic of a periodic tone of 110 Hz.

so the spectrum of $x(t)$ is about what actual frequencies the energy of a signal is.  it doesn't matter whether $x(t)$ is periodic or not (actually it does with the esoteric mathematics, but we won't go there).  and we use the Fourier Transform, not the Fourier Series to visualize it.

So each harmonic in the Fourier Series above is attached to a dirac delta function, $\delta(f)$.

$$\begin{align} X(f) &\triangleq \mathscr{F} \Big\{ x(t) \Big\} \\
 &= \int\limits_{-\infty}^{\infty} x(t) \, e^{-j2\pi f t} \, \mathrm{d}t \\
 &= \mathscr{F} \Big\{ \sum\limits_{k=-\infty}^{\infty} c_k \, e^{j(2\pi k/P)t} \Big\} \\
 &= \sum\limits_{k=-\infty}^{\infty} c_k \,\mathscr{F} \Big\{  e^{j(2\pi k/P)t} \Big\} \\
 &= \sum\limits_{k=-\infty}^{\infty} c_k \, \delta(f-k\tfrac{1}{P}) \\
\end{align}$$

If you express the spectra of your periodic signals in this form, the spectra sum to a single spectrum.


$x[n]$ after sampling of $cos(16\pi t+\phi)$ at 12kHz I'm not sure what the question really means, so this is just guesswork.

I think options 1 and 4 can be ruled out as $w_0<\pi$.

The CTFT of $cos(16\pi t+\phi)$ has two spikes at $16\pi$ and $-16\pi$. I guess that the question means that we sampled this spectrum at $12Hz$, (because 12kHz doesn't make sense to me, might be a typo in the question), i.e. at $24\pi$.

So we get one copy of $-16\pi$ at $-16\pi+24\pi=8\pi$. Similarly we get a copy of $16\pi$ at $-8\pi$.

Now we normalize to the $w$ by diving by $12$. Only the frequencies $-8\pi$ and $8\pi$ get mapped in the range $-\pi$ yo $\pi$ as required by the question. These frequencies get mapped at $-\frac{2}{3}\pi$ and $\frac{2}{3}\pi$. But the problem is that they also get scaled by $12Hz$ during sampling (and I'm ignoring that to get to one of the options). So the $x[n]$ corresponding to $\frac{2}{3}\pi$ and $frac{2}{3}\pi$ is $cos(\frac{2}{3}\pi n+\phi)$

Is this correct?

First of all, as you said the sampling rate is probably $12$ Hz, rather than $12$ kHz, and perhaps they want to demonstrate an aliasing example.

Given a bandlimited continuous-time periodic signal
$$x(t) = \cos(16\pi t + \phi)$$
the samples taken at the rate $F_s = 12$ Hz will be denoted as $x[n]$ and will be obtained by via $x[n] = x(t_n)$ with $t_n = n T_s = n/Fs$ :

$$x[n] = \cos( 16 \pi \frac{1}{12} n + \phi) = \cos(\frac{4 \pi}{3} n + \phi)$$

$x[n]$ is a discrete-time periodic sequence with frequency $w_0 = \frac{4\pi}{3}$, however, in the discrete-time case we are interested in the range of frequencies that fall in $ [-\pi,\pi)$, and therefore this signal's frequency will be cast as:

$$ x[n] = \cos(\frac{4 \pi}{3} n + \phi) = \cos(\frac{2 \pi}{3} n + \phi)$$

which indicates that the original continuous-time frequency of $8$ Hz is lost, and a new aliased frequency of $4$ Hz is attained...


Probability error for antipodal systems and orthogonal systems I know that the error probability for antipodal systems is  Q (√((2Fs)/No)) and for orthogonal systems it is Q (√(Fs/No)), however, can anybody give their derivation? Because I am finding different ones on the internet I am getting confused.

I found pages 74-75 useful but I think their derivations are quite short here:
http://www.sps.ele.tue.nl/members/F.M.J.Willems/TEACHING_files/5JK00/signalenergyorthogonalsignals.pdf

For equally likely
binary signals $s_0(t)$ and $s_1(t)$ of energies $E_0$ and $E_1$ respectively
in an AWGN channel with two-sided power spectral density $\frac{N_0}{2}$, the
error probability is given by
$$P_e = Q\left(\sqrt{\frac{E_0 + E_1 - 2 \langle s_0, s_1 \rangle}{2N_0}}\right)
\tag{1}$$
where
$\displaystyle \langle s_0, s_1 \rangle = \int_{-\infty}^\infty s_0(t)s_1(t)\,\mathrm dt$ is the inner product of the two signals.

For antipodal signals, $E_0 = E_1 = E$ and $\langle s_0, s_1 \rangle = -E$ giving
everyone's favorite answer $Q\left(\sqrt{2E/N_0}\right)$ while for equal
energy orthogonal signals, the inner product is $0$ giving everybody's second
favorite answer $Q\left(\sqrt{E/N_0}\right)$. For some details of the proof
of $(1)$,
see this answer that I wrote
some time ago.  For even more detail,
see
this Lecture Note of mine.


Continuous EQ defined by a curve I'm trying to figure out how to implement a continuous EQ filter.  I have in the past built a parametric EQ and then built a 10 channel EQ using 10 of these parametric filters.  This gave nice results.

However I'm now trying to create a continuous EQ filter.  Is this possible with a parametric EQ?  By selecting the Q factor, perhaps?  If so how would I fit the Q factors to the curve?

Failing that is there another better way to implement continuous EQ?  Could it be implemented from an FFT by extending the magnitude of each r,i bin based on the curve?  Or would this produce rather nasty EQ'ing?

Any thoughts appreciated.

Edit: To answer the comment and add a bit more.

I've been experimenting with trying to define an FFT filter.  Basically I define an arbitrary curve.  The curve defines dB increase or decrease.

I'm currently building my own FFT by looking up the y intercept on the curve for each frequency bin. (Well i'm converting from dB to a multiplicative factor).

Then, later, I'm FFT'ing a block of audio and then I'm multiplying it by the values in the "filter" FFT.  When the multiplicative factor is purely 1 then I get the audio passing through fine.  As I increase the multiplicative factor for a given bin I get a nasty crackling (though I do hear the correct frequency boosted).

Surely there must be a way of doing this in the fourier domain.

Btw, please don't give me matlab functions unless you can explain how they work because I'm trying to implement this on a mobile device so I can't use matlab.

This is a very complicated topic. A good overview on how to do this in MATLAB is here:
https://ccrma.stanford.edu/~jos/vguitar/Fitting_Filters_Matlab.html

An FIR fit can typically done by interpolating up to a full FFT grid and adding either a minimum or a zero phase and then doing an inverse FFT to get the impulse response. You can then try different filter lengths to determine which lengths is required to create an acceptable fit.

IIR fitting is much more tricky. The above example uses MATLAB's invfreqz() which in some cases gives good results and in other only works very poorly. In general the methods are iterative, i.e. there are many design cycles and there are small modifications between design steps that are supposed to "improve" things. These often use "gradient" methods, i.e. you need to evaluate the derivatives of the error function with respect to zeros and poles (or whatever other representation you may choose). They also involve least square error methods. A lot of care has to be taken to correctly formulate the error criteria in terms of frequency spacing, weight functions, relative vs. absolute error etc. The formulation of the error function is a direct representation of what the specific application can tolerate (or not).

EDIT based on OPs question edit:
If you want to do frequency domain filtering, (i.e. FFT->mpy with Filter->inverse FFT) then you need to get all the framing, zero-padding and overlap handling correct. This due to the difference of circular vs. linear convolution and that the FFT "as is" assumes by definition periodic signals. Google "overlap-add" for more info.

You still need to decide on the phase for your filter. Leaving the phase zero creates a non-causal filter that will create lots of wrap-around problem. At the very minimum you need to shift the impulse response so it's causal.


Calculating filter parameters from system response In trying to answer another question, I got stuck applying what I thought would be a valid way to answer it: using the system function. Here's what I did:

First, I put the question

Y[i] = Y[i-1] + ALPHA * ( X[i] + Y[i-1] )

into DF I:

Y[i] = ALPHA * X[i] - (ALPHA-1) * Y[i-1]

This tells me that the system function is:

H(z) = ALPHA / ( 1 + (ALPHA-1) z^-1 )

Substituting:

z=e^(jw)

and

|H(z)| = .5 (or .707)

should allow me to calculate ALPHA, but I only get complex values for ALPHA. What am I doing wrong?

Just math I think. Let's substitute g = ALPHA -1, which makes the algebra a little easier. Then you need to solve
$$\left | \frac{1+g}{1+g\cdot e^{-j\omega }} \right |^{2}=0.5$$
That means
$$\frac{1+g}{1+g\cdot e^{-j\omega }}\cdot \frac{1+g}{1+g\cdot e^{+j\omega }} = 0.5$$
Multiplying this out yields
$$\frac{(1+g)^{2}}{1+2\cdot g\cdot cos(\omega) + g^{2} }= 0.5$$
We can put this into standard quadratic equation form:

$$g^{2}+(4-2\cdot cos(\omega))\cdot g + 1 = 0$$

Solve, plug in the numbers and you get
$$g_{1} = -1.0744, g_{2} = -0.93074$$

The first solution doesn't make sense since it would be a pole outside the unit circle, so the second is correct and we get the expected

$$\alpha = 0.069262$$


How to perform a Stockwell Transform on a 2D matrix I was looking at performing a Stockwell transform on a matrix (or an image as it has two directions - the matrix I am interested in has a spatial dimension and a time dimension, but this should make no difference) instead of an array. Whilst there are a couple of codes that perform a Stockwell transform on an array (see infomation about Stockwell Transform for links to 2 of them and http://www.mathworks.com/matlabcentral/fileexchange/45848-stockwell-transform--s-transform- for one without for loops), I could not find one for a matrix/image.

According to Mansinha et al. (1997), for a matrix $h(x,y)$ the S-transform can be defined as

$ S(x,y,k_x,k_y) = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} H(\alpha + k_x, \beta + k_y) e^{-\frac{2\pi ^2 \alpha^2}{k_x^2}} e^{-\frac{2\pi ^2 \beta^2}{k_y^2}} e^{i2\pi (\alpha x + \beta y)} \mathrm{d}\alpha \mathrm{d}\beta$

for $k_x \neq 0 $ and $k_y \neq 0 $; where $H(k_x,k_y)$ is the 2D Fourier transform of the matrix $h(x,y)$. $x$ and $y$ are spatial variables and $k_x$ and $k_y$ are the corresponding wavenumbers.

The discrete form of the above according to the same paper is given by

$ S(pT_x, qT_y, \frac{n}{NT_x}, \frac{m}{MT_y}) = \sum\limits_{n'=0}^{N-1}\sum\limits_{m'=0}^{M-1} H(\frac{n'+n}{NT_x},\frac{m'+m}{MT_y}) e^{-\frac{2\pi^2 n'^2}{n^2}} e^{-\frac{2\pi^2 m'^2}{m^2}} e^{\frac{i2\pi n'p}{N}} e^{\frac{i2\pi m'q}{M}}$

for $n \neq 0$, $m \neq 0$; where $H$ is the 2D Fourier transform of the matrix $h(pT_x,qT_y)$ and $p=0...,N-1$, $q=0...,M-1$ and $T_x,T_y$ are the sampling intervals in the respective directions.

Therefore (if I understand it correctly) the above process is a inverse Fourier transform to a multiplication in the frequency domain.

Unfortunately, I was unable to code up because the 2D Fourier transform of $h(pT_x,qT_y)$ has the same dimensions as $h$ and thusthere is a mismatch with $H(\frac{n'+n}{NT_x},\frac{m'+m}{MT_y})$.

Any help in coding the above in MATLAB is greatly appreciated. This (http://www.mathworks.com/matlabcentral/fileexchange/45848-stockwell-transform--s-transform-) is a very efficient 1D transform which has helped me understand how it works. The link for the paper is http://www.sciencedirect.com/science/article/pii/S0031920197000472 (if you have access to it)

There is a MATLAB code available for the Stockwell transform here. It also includes a code for the one-dimensional transform but it is less efficient that the Mathworks file exchange code that is mentioned in the question. This is because one uses for loops and the other doesn't.

Using this 2D S-transform function, you can quickly run into memory issues. You can easily modify the code to specify the exact positions of x and y that you need, so that S is not 4 dimensional.


Designing filters vs. editing fft components What is the difference between designing high-complex filters to reduce/remove frequencies in a signal and calculating the fft, reducing/removing those frequencies manually, and reconstructing the signal? Although the latter seems to be easier, I am pretty sure there must be compelling reasons to use filters instead of editing frequencies.
Thanks.

I'll here give another reason (perhaps the most important one) why implementing filters by manually changing frequency bins is usually a bad idea.
Note first that time domain convolution (simple FIR filter design) can also be implemented in the frequency domain and this procedure is known as FFT convolution ( despite its name the math operation here is multiplication and not convolution). The two techniques are equivalent but the latter is much trickier to implement correctly (although it's much faster than its time domain counterpart)
On the other hand, manually changing frequency bins (either by removing them entirely or just by reducing their amplitude),  means that you're discarding the original phase information, which in turn, can lead (and usually does) to inaccurate reconstruction of the modified time domain signal. Multiplying in the frequency domain does correspond to time domain convolution but, unlike true FFT convolution, we don't have a real filter kernel here so the results are less predictible.
Mathematically, perfect reconstruction is possible only when both original magnitude and phase information is retained. But since you're modifying both magnitude and phase, it's quite possible that there is no real signal (in time domain) that corresponds to your modified FFT data (at least not without significant digital artifacts being added to the reconstructed signal). That said, there are ways (e.g iterative procedures such as the Griffin-Lim algorithm ) that can be used to remedy this situation.


Filtering an Audio File in MATLAB I am trying to process an audio file in Matlab by filtering out all frequencies except those within $\pm 25\ Hz$ of $523\ Hz$ (as well as its harmonics up to the Nyquist).

The resulting waveform should look like the green wave displayed below (blue being the original):

However, my resulting waveform is entirely zero:

Here is the simple code I wrote:

[s, Fs] = wavread('file.wav');
n = length(s);
f = linspace(-n/2, n/2, n)/(n*Fs);
S = fftshift(fft(s)/n);

filter = zeros(size(S));
f0 = 523;
while (f0<(Fs/2))
    filter(abs(f-f0)<25)=1;
    filter(abs(f+f0)<25)=1;
    f0 = f0+523;
end

s_filter = real(ifft(fftshift(S.*filter)*n));
s_filter = int16(s_filter);

Tspan = n/Fs;
t = linspace(0,n,n)/n*Tspan;

plot(t,s, t, s_filter);

The problem is that the filter array returns only zeros. But what exactly causes that?

Any help would be greatly appreciated.

First of all, you don't need to shift the spectrum of the signal. Also, your filter is zero everywhere. Your frequency vector contains normalized frequencies, which is fine, but then you are using real frequencies when you build the filter. Either change the frequency vector f to real frequencies or normalize f0 too


Whether Poisson noise will be dominant source of noise in low light conditions? I have one paper saying that that the Poisson noise will be dominant source of noise in low light conditions and the other one says the opposite way(Poisson noise will be dominant except in low light conditions). Which is true?
I believe Poisson noise is due to the random arrival of photons. So it should be more at short exposure time irrespective of light conditions.

The Hasinoff paper says:

  For small photon counts, photon noise is generally dominated by other
  signal-independent sources of noise.

From my own experience with cameras I know that that is true. Even if the lens cap is on, the pixels will have dark current noise that is both spatially patterned and to a lesser extent varies by time. See for example this Nikon D300 dark frame:

The Alter et al. paper only considers the low-light case. They do not take into account any time-varying component of the dark current noise. They do mention an "offset voltage", and take it into account in their noise model. Maybe their sensors are just very good with an almost perfectly constant offset voltage for each sensor element.

The dark current noise is not photon noise, because by definition there are no photons, but it is still Poisson-distributed (see Justin Charles Dunlap, Characterization and Modeling of Nonlinear Dark Current in Digital Imagers, 2014), independent of the signal. Because there are these two types of noise sources with Poisson statistics, I don't think it is wise to call photon noise the Poisson noise like Hasinoff does.


Implementing Convolution in Frequency Domain? Suppose, we have a bitmap image represented as a 2D integer array,
 int [,] image2D; whose FFT is Complex[,] fftImage2D;

Suppose, we have an kernel represented as a 2D integer array,
 int [,] kernel2D; whose FFT is Complex[,] fftKernel2D;

We know that, the convolution (in spatial domain) of image2D and kernel2D would be,

int Rows = image2D.GetLength(0);
int Cols = image2D.GetLength(1);

for(int i=0 ; i<Rows ; i++)
{
    for(int j=0 ; j<Cols ; j++)
    {
        //sweep the kernel2D across image2D
        //...........................
    }
}

We also know that, convolution in frequency domain would be, multiplication between fftImage2D and fftKernel2D.

How can I do this multiplication?

How can I multiply two Complex [,] type 2D arrays of different dimensions? I have understood the theory. My problem is practical implementation. As I described in the question,

Are DFT of the image and DFT of the kernel going to be of different sizes? I guess so. So, how can I multiply them element by element?
In my code, each of the DFTs are represented by 2D Complex numbers. Should, I multiply them according to complex-number's multiplication rule? Probably yes. But, only when their dimensions are same. Right?

IDFT the smaller or both of the DFTs if needed.  Zero pad one or both of the kernel and image to make them the same dimension and size.  Re-DFT as needed, and now you can complex multiply the 2 DFT arrays element-by-element because they will now be the same size.


How to calculate binocular overlap I have two cameras of which I know focal length, principal point, rotation and position. The cameras are in an ideal parallel stereo system.

I want to know with as few calculations as possible how much the images of these cameras overlap (let's say a percentage).

I am trying to follow a guide called Creating correct stereo pairs from any raytracer, but I cannot find how to set aperture and eye separation to find delta.

Can I do it with my current data? Or is there a better way to know how much the camera images overlap?

I'd recommend looking at image stitching algorithms. Even though you don't want to actually combine the images into one, that's still essentially what you're doing.

You might even be able to use an off-the-shelf stitching algorithm to get a rough measure. e.g. If the stitched image is 100 pixels wide, and your two source images are 80 pixels wide, then horizontally you have a 60 pixel overlap.