Playing/Recording acoustic signal in matlab to tell difference between 2 sounds I will be using one central node/latop to play a tone and will have a second node/laptop to record and read the sound and be able to say this node is on the left or right. So trying to figure out the best way to distinguish between the two tones. I want to have a tone set for left and one set for the right. I want the sound to have 400 Samples at 44100 for 10ms between 18-22freq just outside the human hearing range so that the sound is not annoying etc. I want the sound to be different slightly but be able to tell the diff between them for the system to say the node that recorded and processed the sound is on the left of right. How would I go about doing this in matlab because I don't dont really code in this language much and only time I did was for image processing just doing this so I can graduate.

I understand that you want to tell if your second laptop is placed on the left or right (or maybe on the axis) of main node?

I suggest you then to play the two-tone stereo signal, i.e.: $15\;kHz$ from the left, and $16\;kHz$ from the right speaker at the same time. Assuming that you are going to use sampling frequency of $f_s=41000 \;  \text{Hz}$ and taking the average distance between two speakers in a laptop to be $25\;cm$, you can calculate that the maximum delay between two signals, if the second node is on the right, will be approximately: $\dfrac{0.25}{343}=0.7289\;msec$. This corresponds to $32$ samples at a given sampling frequency. It is therefore a hard task, but still possible. Then after finding the delay between two signals, from its sign you are able to tell which speaker is further away, ergo you are on the opposite side.

In MATLAB you just need to generate proper signals, given sampling frequency and play them back/record with: MATLAB - Record and Play Audio.

Few problems can arise:

You want to use very high frequencies, and this can be problematic as laptop speakers are usually low-grade gear and they possibly might have some problems- same goes to microphone. Also remember to turn any hardware noise cancellation, spatial enhancements, etc.
Detection of you signal onset could be problematic as the resolution of DFT will be to low. You can also use some kind of correlation detection, but this is a problem for later on.


Effect of different downsampling rate on EEG signal I have a raw EEG signal of length 2,50,000 samples. I band passed it from 0.05-10Hz (sampling frequency=500). Then I down sampled it at 20Hz. Now question is: What will be effect if I vary down sampling from 20Hz to say 5Hz or 60Hz? And why 20Hz is appropriate?

Given that you have removed all signal content above 10 Hz, your 20 Hz sample rate will be identical to the minimum sample rate based on the Nyquist sampling criterion.

If you instead use 5 Hz, you can only represent a bandwidth of 2.5 Hz and, e.g., a signal appearing to be at 2 Hz might have been a 3, 7 or 8 Hz signal as well (meaning that you can not distinguish between a signal at any of those frequencies after the downsampling). This is caused by aliasing.

Using 60 Hz would be clearly feasible but will consume more storage of the signal. The possible benefit might be that later processing may be simpler as the signal does not cover the whole signal band. However, this depends on which processing is required later. Also, while the processing in itself may be simpler, there are more samples to operate on so the total complexity may not be reduced.


Techniques for variance stabilizing? Filtering Poisson noise? This question is solely theory related so there isn't any sort of practical implementation I'm doing.
I'd like to know what are the different techniques for variance stabilizing(not necessarily used for Poisson noise), and if there's somewhere I can read some detailed explanations of those.
Also I am interested in methods used for removing Poisson noise (the one I found is Ascombe transform, which is technique for variance stabilizing at the same time), again the same goes here, if there's some detailed explanations of these it'd be really great.

Here is the link explaining variance stabilizing transform.

Some papers on Poisson noise removal

Efficient gradient projection methods for edge-preserving removal of Poisson noise
Wavelets,Ridgelets and Curvelets for Poisson Noise removal(Anscombe transform is used here)
Unsupervised Non Local Means for Poisson Noise

Similar question has been asked before: How to remove Poisson and Gaussian Noise?


Understanding the lowpass filtering of a digital signal I do not have any background in Electrical Engineering and I have recently started a project which involves signals captured from sensors.

Lowpass filtering: The way I understand it is: The values below the CUTOFF are allowed to pass through, above the CUTOFF are simply filtered.

So far so good. Now, how does this work in real world?

The data I have is captured 100 times per second at equidistant intervals, hence my sampling frequency is 100 Hz. This is also understandable.

But now, my CUTOFF is 20 Hz. What does it mean? Because all I have is some digits. None of them is below 20. When I read about filtering it makes me think of CUTOFF as some value from the range of signal itself. I have given multiple explanations to myself.

Somehow, from the signal, find out the frequency of each data point and if that frequency is below 20 Hz, don't do anything. Otherwise put 20 Hz in the output. But I have no idea how to find out frequency of single data point. After all, its just a number, right?
It has to work on multiple values at a time. That is, from one data point to the other, if there is some above average peak, just smooth it out. And this 20 Hz is essentially a limit of how sharp the peak can be.

Please help me understand in the simplest terms. I asked the same question to my guide and he explained me using the terminologies in Electrical Engineering and some more alienate stuff, which was ways and bounds beyond my brain.

Any signal can be constructed from the sum of sinusoids of different frequencies (to within an arbitrarily small error).  This reconstruction is unique and can be calculated with the Fourier transform.

When we talk about the various frequencies that comprise a signal, we are talking about these sinusoids of various frequencies.  In fact, in signal processing, we are used to thinking about the same signal in the time domain -- as a function that gives a magnitude for every instant in time, or in the frequency domain -- as a function that gives a magnitude and phase for every frequency.

When you apply a low-pass filter, you are modifying the frequency domain view of the signal so that all frequencies higherthan the cutoff have their magnitudes set to pretty much zero.

Another way of thinking about it is, you take all the sinusoids that make up the original signal, and throw away the ones with frequency higher than your cutoff.

The secret to how this actually works is called the "convolution theorem": https://en.wikipedia.org/wiki/Convolution_theorem

The convolution theorem shows that convolving two signals together is the same as multiplying their frequency-domain representations together.  A "low pass" filter is just a signal that has a constant magnitude for frequencies below the cutoff, and near zero magnitude for frequencies above the cutoff.

When the filter is convolved with the input signal, the input signal's high-frequency components are multiplied by zero, effectively removing them.


Are discrete systems defined by LCCDE always LTI? Suppose a discrete-time system is defined by linear constant-coefficient difference equation

$$\sum_{k=0}^{N} a_k y[n-k] = \sum_{k=0}^{M} b_k x[n-k]$$

where at least two different coefficients $a_i,a_j$ are nonzero (so the equation is recursive). Isn't it true that such system is LTI and causal regardless of any auxiliary condition? After all, it is causal because it does not depend on any future inputs (RHS of the equation has only present or past inputs). And it is LTI because: if $x_1[n] \to y_1 [n]$ and $x_2[n] \to y_2 [n]$ the above difference equation is valid for $x_1,y_1$ and $x_2,y_2$ and by adding them together we can show that

$$\sum_{k=0}^{N} a_k (y_1[n-k]+y_2[n-k]) = \sum_{k=0}^{M} b_k (x_1[n-k]+x_2[n-k])$$

which means $x_1[n]+x_2[n] \to y_1[n]+y_2[n]$.

(similarly we can prove the homogeneity and the time-invariance property).

However it doesn't seem to be true according to Signals & Systems by Oppenheim A.:

  As in the continuous-time case, such equation does not completely specify the output in terms of the input. To do this, we must also specify some auxiliary conditions.

Up to this point I agree. But then they write:

  [..] we will focus for the most part of the condition of initial rest - i.e., if $x[n] = 0$ for $n < n_0$, then $y[n]=0$ for $n<n_0$ as well. With initial rest, the system described by the [first] equation is LTI and causal

In other words, the author seems to imply that without the initial rest condition, there's no guarantee that the system is LTI and causal (why contradicts my previous statements). Where's my mistake?

What you have is a so-called behavior which truncates the signal space $\mathcal{X}\oplus\mathcal{Y}$ of all possible signal pairs $(x, y)$ to the ones that satisfies this equality.

This does not specify an input-output mapping and hence a dependence and hence causality does not apply here. They simultaneously satisfy this. However regardless of causality, it is LTI since you can take arbitrarily and superposition applies.

Just like you can integrate any differential equation sufficiently many times and end up with a integral equation, you can also forward time shift these equations then if one of them depends on future inputs and that is your choice of dependent variable then you can speculate about causality.

This is actually the tool to allow for nondifferentiable input signals and other discontinuities and the results are investigated under the name of weak solutions


HRTF implementation in matlab I have set of HRIR( Head related impulse response) downloaded. I want to convolve it with an audio signal to generate 3D signal in Matlab.

I am new to this field. Is there any tutorials or can any one help me how to start

Each HRIR is meant to represent the acoustic impulse response from an audio source at a specified angle (relative to the head) to the ear canal.  There is no such thing as a 3D signal, but if you want to simulate audio from a given azimuth and elevation angle, just convolve your audio source with the corresponding HRIR as follows:

left_channel = conv(audio_signal,hrir_left_az_elev)
right_channel = conv(audio_signal,hrir_right_az_elev)

In this case, I am using hrir_left_az_elev to signify the HRIR that corresponds to the left-ear impulse response to a sound event at azimuth angle az and elevation angle elev.  How you actually code this up will depend completely upon the format of your HRIR data.


Decorrelating Stationary Colored Gaussian Noise -- Effect On The Desired Signal So given stationary colored gaussian noise $\mathbf{n}$, I know that I can decorrelate it by first finding it's autocorrelation $R_{nn}$ and performing $R^{-\frac{1}{2}}_{nn} \mathbf{n}$.

In practice of course, I need to estimate $R_{nn}$, which I can do via averaging $\mathbf{r}_{nn}$ or by fitting it to an AR model.

So, suppose I have a good estimate $\hat{R}_{nn}$, for a finite data record of length $N$. Such that the resulting spectrum of the noise frame is white enough, according to some metric.

Now, suppose I have a signal $\mathbf{x} = \mathbf{s} + \mathbf{n}$. I want to decorrelate the noise, such that after the decorrelation, I get a signal in white noise. From what I have read, the way to do this is to also apply $R^{-\frac{1}{2}}_{nn} \mathbf{x}$.

However, won't this distort the desired signal $\mathbf{s}$? Can anyone point me to some references regarding this?

Also: How else can I decorrelate colored noise in order to get my signal embedded in white noise?

Yes, the "Whitening" is basically filtering the signal using an LTI System.

The Result would be the signal $ \mathbf{s} $ filtered by the system which whitens the noise.

In the framework of Matched Filter this is OK since you will be able to filter the Reference Signal as well and hence have the same performance.

Yet for other usages you must take under consideration the effect of whitening system.


Sending complex signal over a single wire The G.hn protocol uses OFDM to transmit data. And it sounds like they are doing it on a single wire. OFDM symbols are a complex number with I and Q components. How are they able to send two separate data streams on a single wire?

If they are not doing it on a single wire, then my question is really how do wired protocols send complex I/Q data over a single wire? I'm looking for a fairly concrete answer (does not have to be specific to G.hn) as I'm interested in implementing something like this myself along a single channel.

Note that signals sent over wires (and over the air, and over any medium) are always real. What quadrature means is that, on a passband channel (wired or otherwise), you can transmit two signals at the same time. The first signal, which we'll call $s_I(t)$, is mixed with a carrier $c_I(t)=\cos(2 \pi f_c t + \phi)$; the second signal, $s_Q(t)$, is mixed with $c_Q(t)=\cos(2 \pi f_c t + \phi + \pi/2)$, where $f_c$ is the center frequency of the channel.

Note that it is important to use a passband channel. Quadrature communication is impossible on the baseband channel. Maybe this is what confused you, since one tends to associate wires with baseband. However, as long as you don't exceed the bandwidth of the wire, it's perfectly possible to transmit modulated signals over them. In fact, in telephony, the very first trunk links used single sideband modulation to multiplex several phone calls over a single wire -- this was before PCM and digital telephony were even invented.

The transmitter would then calculate the following signal: $$s(t) = s_I(t)c_I(t) - s_Q(t)c_Q(t),$$ where the negative sign is just a convention. Assuming the carrier phase $\phi=0$ for simplicity, this signal can also be written as $$s(t)=\text{Re}[(s_I(t)+js_Q(t))e^{j 2 \pi f_c t})].$$ This expression is sometimes useful in a digital implementation, since you only have to keep track of two complex signals instead of four real signals.

The channel adds noise, so the received signal is $r(t) = s(t) + n(t)$. The noise is usually modeled as Gaussian, white noise with PSD equal to $N_0/2$.

In the receiver, we have (ommitting a lot of algebra readily found in textbooks and on other answers in this site): $$\hat{s}_I(t) = s_I(t) + n_I(t) = \text{LPF}[r(t) c_I(t)]$$ and $$\hat{s}_Q(t) = s_Q(t) + n_Q(t) = \text{LPF}[r(t) c_Q(t)],$$ which means that the receiver can recover both transmitted signals (albeit with noise).

How to actually get all this to work in an actual system requires solving other problems that have not been addressed here and that can take a whole textbook to explain. I recommend that you start with this free textbook: http://sethares.engr.wisc.edu/telebreak.html (link at the end of the page).


How to estimate octave and size for visual features positioned at Harris corners I am currently working with and comparing the performance of several feature detectors provided by OpenCV as basis for visual feature matching.

I am using SIFT descriptors. I have accomplished satisfactory matching (after rejecting bad matches) when detecting MSER and DoG (SIFT) features.

Currently, I am testing my code with GFTT (Good Features to Track - Harris corners) to get a comparison, and also because in the final aplication, a set of GFTT features will be available from visual feature tracking process.

I am using cv::FeatureDetector::detect(...) which provides me with a std::vector<cv::KeyPoint> filled with detected features / keypoints / regions of interest. The structure cv::KeyPoint contains basic information about feature's location, as well as information about size and octave in which the keypoint has been detected.

My first results with GFTT were terrible until I compared the typical size and octave parameters in different types of features:

MSER sets the size (between 10 and 40px) and leaves the octave to 0
DoG (SIFT) sets both the size and the octave (size/octave ratio between 20 and 40)
GFTT the parameters are always: size = 3, octave = 0

I presume that is because the primary purpose of GFTT features was not to be used in matching but only in tracking. This explains the low quality of matching results, since the descriptors extracted from such tiny features stop being discriminatory and invariant to many things, including small, 1-pixel shifts.

If I manually set the size of GFTT to 10 - 12, I get good results, very similar to when using MSER or DoG (SIFT).

My question is: is there a better way to determine how much to increase the size (and/or octave) than just-go-with-10-see-if-it-works? I want to avoid hardcoding the size increase if possible and determine it programmatically, but hardcoding is okay as long as I have some solid arguments backing up my choices of the new size / size increase / size estimation algorithm.

I'm not sure there is actually a good response to your precise question : the scale-space thing of SIFT and SURF was actually developed to automatically estimate the "good" relevant neighborhood size around a corner-like keypoint (which is what good features to track are).

Now, more positive answers would be:

build a database of keypoints and good matches (e.g. using square calibration patterns) and create an automated performance assessment on this database to find the correct size. This task can actually really be automated (see the work by Mikolajczyk and Schmid about point matching evaluation)
embed your good features in an image pyramid to also have some kind of scale associated with them. You can look for references on multi-scale FAST and Harris interest points, that do something very similar to this procedure.

To heuristically find the maximum block size, you can compute estimates of your images with a box blur (which is more or less what the blockSize operator does) and see when the corner disappears. Note however that more blur takes the corner away from its true location.

If you're really looking for some quick-and-dirty fix, try sizes between 5x5 and 11x11 (typical sizes used in stereo block matching). If you're looking for an intellectually satisfying criterion, then try to maximize the likelihood of good matching of two feature points under your noise level.


difference equation block diagram delay or advance confusion When does the unit delay, become unit advance essentially.
Goal is to deduce difference equation from block diagram.
Green text is notes about what the value should be in that node of the diagram.
It appears that the difference equation should be equal to
\begin{align}
y[n+2]&=x[n]&-&0.9y[n]\\
x[n]&= 0.9y[n]&+&y[n+2]
\end{align}

Green text denotes my personal notes.

It depends from which side of the delay block you look at. Below is an example one sample delay from left to right in the arrow direction :

$$ x[n] \longrightarrow \boxed{ z^{-1} }  \longrightarrow y[n] $$

Then either $y[n] = x[n-1]$ or $x[n] = y[n+1]$ are both correct and equivalent expressions.

The first one is written when you trace the delay block in the arrow direction. And you get a delay effect from input x to output y. This is the normal operation direction.

The latter expression is obtained when you trace the delay block in the opposite direction of the arrow. Which gives you an advance effect from signal y to signal x.


How can I detect if a filter is FIR or IIR according to impulse response graph? I have the following two graphs but I don't know how to recognize from the graph if they are IIR or FIR! Can you please guide me through?

P.S.: Does MATLAB have a function to detect if an impulse response graph corresponds to IIR or FIR?

If a filter is FIR then the impulse response is finite and goes to zero and stays there. This is like the second filter in your post. If it is IIR the impule response extends towards infinity like the first image in your post


How can I find the amplitude of the 20 Hz frequency component of a signal faster? I have a time varying discrete signal x(t) sampled at 5ms rate.
I need to find the amplitude of the 20 Hz frequency component for each sample.
What I do today is that for each sample I do a Short-time Fourier Transform of the data in a window around that sample. Then I read the amplitude for the 20 Hz frequency component off the  amplitude spectrum.
This seems to work fine, but it is slow.

How can I speed up this process?
For instance I am only interested in the 20 Hz component. Is there some way I can avoid finding the complete amplitude spectrum (that is faster than FFT)?

As pichenettes already pointed out in a comment above, one approach is to use the Goertzel algorithm. However, it is a somewhat common misconception that in order to deduce the amplitude of various frequency bands in a signal, one must use Fourier-transform-based techniques. There are other methods that may be more appropriate to your application.

In your case, for example, you would do equally well by just creating any bandpass filter centered at 20 Hz and watching the amplitude of its output over time. The DFT (and by extension the Goertzel algorithm) can be viewed as just a bank of equally-spaced, critically-sampled bandpass filters. If you only need one output, then a DFT is overkill, and the Goertzel algorithm can be limiting, as you don't have much control over its frequency response about 20 Hz.


Amplitude Shift Keying coding I'm a computer science student and I took a network course this semester which assumes some signal processing backgrounds which I have not. Last time We talked about ASK.

ASK is Amplitude Shift Keying,

it uses two different amplitudes to describes $0$ or $1$. We can use several levels to represents more symbols.

$$A\sin(2\pi fnt+\phi)$$

How many harmonics does a signal with an ASK coding has?

If I have such a message to broadcast: $01101101$

How do we code in ASK this series of bits when the frequency used is: $f=1kHz$? And when the bit-time is $T_b=1ms$?

And what about varying its amplitude? If we have 4 levels for instance?

Regarding the harmonics: ASK has an infinite number of them. The ASK signal can be thought of as a square wave multiplied by a sinusoid. In frequency, this corresponds to the spectrum of the square wave shifted to the frequency of the sinusoid. Since the spectrum of the square wave is infinite, the spectrum of the ASK signal is infinite too.

In practice, the ASK signal is filtered to force it to fit in a specific bandwidth. This has the effect of distorting the signal, but if done carefully the distortion is not severe. In any case, ASK is only used in applications where bandwidth is not a constraint.

I recommend the following procedure to understand how to generate an ASK signal. First you need to map "bits" to "voltages" (also called "amplitudes"). Bits are abstract quantities that do not have physical substance, while voltages are physical. In ASK you typically map a bit 0 to 0 volts, and a bit 1 to $A$ volts. Now you can generate a sequence of voltages that corresponds to your bits:

$$ \text{bits}=\lbrace 0,1,1,0,1,1,0,1 \rbrace$$

$$ \text{amplitudes} = \lbrace 0, A, A, 0, A, A, 0, A \rbrace $$

Now you need to define a pulse $p(t)$ of duration $T_b$:

$$ p(t) = \begin{cases}1,\quad0\leq t\leq T_b,\\0,\quad\text{otherwise}\end{cases} $$

(You will eventually want to control the energy of the pulse by changing its amplitude).

The next step is to define the "baseband" ASK (that is, with no carrier) signal:

$$ s_{BB}(t) = \sum_{k=0}^7 a_k p(t-kTb) $$ where $a_k$ is the $k$-th amplitude (starting from 0), and $p(t-kT_b)$ is the pulse shifted $kT_b$ seconds "to the right". Take a moment to study this equation and you'll see that it is a train of pulses of amplitude 0 for a bit 0 and amplitude A for a bit 1.

Now you're ready to create the ASK signal, which is:

$$ s(t) = s_{BB}(t)\sin(2\pi f_c t) $$

To transmit more bits per pulse, all you need to is to increase the range of values of the amplitudes.

I know that this is different from the way most undergrad textbooks teach this subject, but I still think it makes more sense to think about it this way.


Quantization Image using MATLAB I'm trying to quantize an image 8 bits to 4 or 2-bits uniformly. I searched internet, interestingly I could not find what I want exactly. Then I wrote an simple code for it myself. I'm curious about whether there is a build-in function in MATLAB which convert 8-bit image to 4-bits uniformly . My results using the methods from internet is not good. Am I doing something wrong when I use these methods? Thank you in advance.

Method from internet giving strange results

reducedImage = uint8((single(monalisa)/256)*2^4);

Method writen by me. I'm putting code to clarify what I want to do.

monalisa_2= monalisa;
figure
monalisa_2(monalisa<=63)               = 32  ; %% 00
monalisa_2(monalisa<=127&monalisa>63)  = 85  ; %% 01
monalisa_2(monalisa<=191&monalisa>127) = 159 ; %% 10
monalisa_2(monalisa<=255&monalisa>191) = 223 ; %% 11
imshow(monalisa_2)

I point you to this https://stackoverflow.com/questions/12723699/changing-image-bit-depth-using-matlab it can be done. This particular post is only for png, but if I'm not mistaken all image types have the depth parameter so they can likely all be changed in a similar manner

EDIT

For your particular problem, since your images are already grayscale you would do something like

%change this as needed
desired_bit_depth = 2;

my_pixel_depth = 2 ^ desired_bit_depth;
%converts grayscale to an indexed image at the appropriate depth
[ind_im, reduced_colormap ] = gray2ind(my_gray_image, my_pixel_depth)

%displays your new image
figure()
subplot(1,2,1);imshow(my_gray_image); title('original')
subplot(1,2,2);imshow(ind_im,reduced_colormap); title(sprintf('reduced to %d bits',desired_bit_depth))

% save indexed png
imwrite(ind_im, reduced_colormap, 'test.png', 'bitdepth', desired_bit_depth);

EDIT #2: While reducing the "bits" to display the image, the "colormap" is made by taking the average of each range
I kept the original answer because it still works and may suit someone elses needs in the future.

For this solution, rather than using evenly spaced values for our "colormap' we actually use the gray value that is the average over the specified range

ex

to convert to 1 bit image we have two ranges 0,128 and 128,256

replace all values between 0,128 with the average value between 0,128

replace all values between 128,256 with the average value between 128,256

function result_im = ChangeBitDepthGrayImage(gray_im, desired_bit_depth)
    if (desired_bit_depth < 1)
        disp('converting to binary(1 bit) image');
        desired_bit_depth = 1;
    end

    if (desired_bit_depth > 8)
        disp('converting to 8 bit image');
        desired_bit_depth = 8;
    end

    %assuming we start with 8 bit image 256 levels
    num_levels = 2 ^ desired_bit_depth;

    %figures out how big each range should be, we use +1 because if we
    %divide the data into N levels, there should be N+1 boundaries
    limits = linspace(0,256,num_levels + 1);

    result_im = uint8(zeros(size(gray_im)));

    for i = 1:num_levels
        lower_lim = limits(i);
        upper_lim = limits(i+1);

        %creates a binary mask of values between the limits, the output is
        %0 or 1, but we need to make it uint8 for the next step
        temp_mask = uint8((gray_im >= lower_lim) & (gray_im < upper_lim));

        %multiplies image by mask, this isolates only pixels in the given
        %range
        image_only_in_range = temp_mask .* gray_im;

        %finds the mean of that small part of the image. this weird notation is
        %taking the average of nonzero elements
        avg_val_for_range = round(mean(image_only_in_range(image_only_in_range~=0)));

        %replaces all pixels in that range with the average val
        result_im = result_im +(avg_val_for_range * temp_mask);
    end

    %i just picked some random figure
    figure(32)
    subplot(1,2,1);imshow(gray_im);title('original image');
    subplot(1,2,2);imshow(result_im);title(sprintf('modified to %dbit image',desired_bit_depth));
end


Quantifying goodness of amplitude estimation Assume you have a single sinusoid in bandlimited Gaussian noise with unknown amplitude $A$, known frequency $f_0$, and known noise spectral density $S(f)$ in $\frac{\mathrm{units}^2}{\mathrm{Hz}}$:

$$x(t) = A\sin(2\pi f_0t) + n(t)$$

The signal is sampled for a known finite duration $T$ such that the frequency component of the sinusoid alone would have finite magnitude $$\bigg\lvert\mathscr{F}\left\{A\sin(2\pi f_0t)\cdot \mathrm{rect}\left(\frac{t}{T}\right)\right\}(f_0)\bigg\rvert = AT$$.

How would one then quantify the 'goodness' of amplitude estimation $\hat{A} = \frac{|X(f_0)|}{T}$?  I assume the SNR would be a distribution with the following shape?  (Since the signal amplitude would be divided by an integral of zero-mean noise, which is also zero-mean)

Edit: Actually, using the definition of SNR as $\frac{E_s}{\sigma^2}$, SNR is constant since we would know both $E_s (=\frac{A^2T}{2})$ and $\sigma^2 (= S(f) \cdot BW)$.  However, now I'm a bit confused.  This would make it seem that reducing the BW to an infinitesimal amount would make the SNR very high.  On the other hand, I thought the fourier transform was essentially applying a very selective BP filter at each frequency but still results in noise corrupting the frequency magnitude response because there is also a finite amount of energy per bandwidth in the 'smeared' (time-windowed) sinusoid.

Edit2: I believe the main sub-problem is Frequency magnitude distribution of noise

Your understanding of the FT as several extremely narrow bandpass filters is fine. However, you mix up the terms SNR and noise variance. SNR is defined by

$$SNR=\frac{E(s(t)^2)}{E(n(t)^2)}$$

where $E[n(t)^2]=\int_\mathbb{R}S(f)df$ due to the Wiener-Khintschin Theorem ($S(f)$ is the PSD of the noise). So, if you have a measured signal $x(t)=s(t)+n(t)$ with $s(t)=A\sin(2\pi f_0t)$ and $n(t)$ is bandlimited AWGN with bandwidth $B$ and density $N_0$, the SNR decreases when the noise bandwidth increases. So, in order to have a good SNR, you should keep the bandwidth as small as possible (as narrow, as your $s(t)$ permits).

However, your problem is to estimate the amplitude $A$ of $s(t)=A\sin(2\pi f_0t)$ where $f_0$ is known to you. Optimally, you take the Fourier transform of your signal and see what is the value of $X(f)$ at the frequency $f_0$. What will you see there? Does this measurement depend on SNR? Does this measurement actually depend on the noise bandwidth? No, the measurement $X(f_0)$ is independent of the noise bandwidth (as long as the noise has contribution at $f_0$). So, by performing bandpass filtering, increasing/decreasing bandwidth you will not get any better estimates.

The only thing which you can do is to increase the duration of the measurement window. If you have a finite length measurement of duration $T$, your frequency domain essentially has resolution $1/T$. Hence, on each frequency bin (assuming the DFT) the noise variance is $N_0/T$ (assuming AWGN). So, the longer the measurement, the better the estimate.

You can also understand this from a different point: What you are interested in, is the value

$$\rho=\frac{\int_{f_0-\Delta_f}^{f_0+\Delta_f}S(f)df}{\int_{f_0-\Delta_f}^{f_0+\Delta_f}N(f)df}$$

where $2\Delta_f$ describes the frequency domain resolution, $S(f)$ is the signal spectrum and $N(f)$ is the noise spectrum. The numerator is always the same (independent of $\Delta_f$) since $S(f)$ is just a (weighted) Dirac at $f_0$. But, the denominator depends on $\Delta_f$.


Why low-pass FIR Equiripple filter is designed with 3 dB of ripple in the pass-band and at least 60 dB of attenuation in the stop-band? Please tell me, why low-pass FIR Equiripple filters are designed with 3 dB of ripple in the pass-band and at least 60 dB of attenuation in the stop-band
and not other values ?
I want to know are these values simply a convention?

Stop band attenuation is a good thing right? After all it's the primary reason why you want a frequency selective filter: stop the unwanted frequencies. So a 60dB attenuation in the stop band of a filter is considered to be the norm in many typical applications. But not always so, as the comments indicate; you must always refer to your application requirements for determining how much attenuation is enough (or necessary) for your case.

That being said, there exist tradeoffs in filter design, just like in many other design procedures. And the most fundamental of those tradeoff parameters, for digital FIR filter design, are namely the passband ripple, stopband attenuation, transition band width and the resulting filter's length.

Among these parameters transition band width and stop band attenuation affects the selectivity quality of your filter. But a short transition band with a large stop band attenuation means a (very) high order (long impulse response) FIR filter, which creates opposition due to computational burden.

To overcome this computational burden, due to narrow transition band, once can consider relaxing the filter's definition in, for example, the pass-band region, so that instead of a monotonic or an un-equal ripple filter, you distribute the ripples equally by allowing as large as possible but fixed amplitude ripples in every where through the passband, so that the FIR filter length will be kept short enough.

3dB passband ripple can be considered a reasonable choice for some applications in which pass band distortions due to those ripples may be tolerable.


Calculating values of frequency bins in Python I have to use FFT to determine the period of waves inside a signal. After applying FFT on a window of 10000 point from a signal, I get something like this:

What I don't understand is that FFT is supposed to return frequencies, but if the input is a longer signal with the same frequencies, the values of frequencies returned by FFT will change.

So for an array of N length, the result of the FFT will always be N/2 (after removing the symmetric part), how do I interpret these return values to get the period of the major frequency?

I use the fft function provided by scipy in python.

Edit:

Some answers pointed out the sampling frequency. I don't understand what the number of samples per second has to do with the size of the periodic pattern, the FFT returns frequencies right? And then for a specified frequency f, I can do t=1/f and then t will be something like 300 points for example. That means we have a repeating pattern every 300 point, I'm not sure about this. To summaries I have a periodic signal and I need to identify the period T0

This is simply how Discrete Fourier Transform (i.e. Fourier Transform theory applied on sampled signal) works.
You get an output of length N if your input has length N, and after removal of symmetric part, what you get are $\frac{N}{2}$ points that span frequencies 0 (DC component) to Nyquist frequency ($\frac{F_s}{2}$).

For this same reason, the more points your signal have (let's say for a fixed sampling rate $F_s$) the more resolved your FFT output is since the spacing between each frequency bin in the output is:
$$
\frac{F_s/2}{N/2} = \frac{F_s}{N}
$$

UPDATE:

Regardless of the sampling frequency the FFT returns $N$ data points for an input with $N$ samples. What we say is that all the points of the output spread the circle: $[0;2\pi]$ in radian/sample or sometimes denoted $[0;1]$ in cycles/sample (called the normalized frequency). All these bins correspond to how many cycles, or oscillations you might have between each samples. You can measure at best half a cycle between each sample (for real signals). This is usually converted back into frequency simply by knowing the sampling rate.

So indeed if your signal is longer, you will have a longer output. Now what I was pointing out is that if your signal is longer but recorded with same $F_s$, i.e. it actually last longer in time, then you will resolve smaller frequency (since you have more time in total, this gives finer information in the frequency domain). That is, the bins are closer to each other in frequency domain ($\Delta f = F_s/N = 1/T$ get smaller). So basically, the longer you record, for whatever sampling rate, the finer resolution you get in the frequency domain.

In summary, if you do not know $F_s$ the bins are located (after removal of symmetrical part):

span $[0,0.5]$, and spaced by $1/N$, and the unit is cycle/sample, OR
span $[0,\pi]$, spaced by $2\pi/N$, unit is radian/sample

Then if you know $F_s$, simply transform back into Hertz, knowing that it must span $[0, F_s/2]$.


Help me understand how Pulse-Amplitude Modulation (PAM) works Reading on the Internet i learned that PAM is basically a sampling technique which converts an analog signal into a discrete signal. It modulates a train of pulses by varying the amplitudes of the individual pulses in a regularly timed sequence based on the analog wave. (Multiplying the analog signal with a train of pulses)

In my mind i have the idea of PAM being similar to AM but instead of an analog carrier we have a train of pulses.

On the other hand in some textbook i read that the way PAM works is by grouping bits from a stream into symbols and each symbol corresponds to a pulse with a certain amplitude.

This seems to me as a completely different process than the first definition.
So what is true? and if both are true when is the second way to do PAM used?

If you define pulse amplitude modulation (PAM) as a method where a train of pulses is generated with the information coded by the amplitude of these pulses then both methods comply with that definition.

In analog PAM, the amplitudes are defined by samples of the analog source signal, whereas in digital PAM, the amplitudes are discrete and are defined by the bit pattern and some fixed mapping between bit patterns and discrete amplitude values.

In both cases, the PAM signal is given by

$$s(t)=\sum_kA_kp(t-kT)\tag{1}$$

where $A_k$ are the amplitudes representing the analog or digital source signal, $p(t)$ is the pulse, and $T$ is the symbol rate.

The signal $(1)$ is a baseband PAM signal, which can be modulated by a carrier to obtain passband PAM. Using complex symbols $A_k$ on a rectangular grid in the complex plane and two orthogonal carriers gives you quadrature amplitude modulation (QAM). Take a look at this answer for more information on passband PAM.


Synchronization of two audio tracks via Mel-frequency cepstral coefficients (MFCCs) I am trying to solve the (common) problem of synchronizing (i.e. measuring the delay between) two audio tracks recorded from different sources during a conference talk (about 1h of audio). One track is recorded with the not-so-good video-camera internal mic (but in sync with video), and the other one recorded directly from the mixer output (good audio quality). The final goal is obviously to superimpose the good audio from the mixer to the video track. By the way, I am trying to implement this on my own in python.

I have already tried to calculate the cross-correlation between selected segments (ranging from 2min to 20min) of the two (normalized) raw audio files, but it appears that the result is very noisy and extremely dependent on the choice of the segment length and position, making it difficult to search for well defined correlation peaks.

I have stumbled into the approach proposed here by @pichenettes and I found it really interesting:

  Extract a sequence of MFCC vectors from both signals. This way you'll
  have something of lower dimensionality and a bit more robust to noise
  and differences in transducers.
  Optionally normalize the MFCC ("by
  design" coefficient 0 has more variance than coefficient 1 and so
  on...)... Ideally you'd want each row of the MFCC matrix to have
  variance 1.
  Compute the cross-correlation of two matrices you get,
  along the time axis (that is to say, compute the 2D cross-correlation
  and just keep the data for the time axis).

However, I have trouble dealing with the result of this correlation. In my understanding, I would then search where the correlation is maximum in terms of frame index, and then calculate the delay from such index. The point is that after calculating the correlation, I am left to deal with this quantity for 13 cepstral coefficients, which also turns out to be strongly dependent on the coefficient itself. My question is: how can I deal with this array of 13 cross-correlations to calculate the track delay? Am I thinking this in the wrong way?

It is probably useful to mention that the delay between the tracks (measured by ear) is about 18-19s. The sync resolution I need would be something slightly below the 25fps of the original video, e.g. 10ms is definitely ok.

Thank you!

I think, you do not need to follow complicated steps of MFCC. I recommend you to preprocess the signal. So you can remain the portion between 300~4000Hz. (1. FFT, 2. assign zeros for the frequencies for smaller than 300 and larger than 4000Hz. 3. IFFT) And then you can do cross-correlation stuffs like what you did. Hope it helps.
You can understand the idea of MFCC from here.

Then you will understand why I don't recommend you MFCC but to remain 300~4000Hz.


Why do I have to use variance and not mean value to calculate noise? Why, when I have to calculate some types of noise (for example shot noise and others which are in the denominator of SNR ratio) or attenuation from turbulence, I have to use variance of the of the density probability and not the mean value? In all of the formulas they use  variance (the second moment) and this is not clear for me.

Typically the first moment (mean) is zero and therefore contains (relatively) no information. The most obvious thing then, to distinguish different levels of noise, is to look at the variance.

And if the variance happens to be zero, then surely there is no noise at all.


Finding the maximum phase and the maximum frequency deviation of an angle-modulated signal I have an angle-modulated signal:
$$
s(t)=10\left[\cos\left(10^8\pi · t + 5\sin\left(2\pi\cdot10^3·t\right)\right)\right]
$$
How would I go about finding the maximum phase and the maximum frequency deviation?

I think this is treating the $10^8\pi\cdot t$ as the carrier frequency and the $5 sin(2\pi\cdot 10^3 t)$ as the modulation.  So the maximum phase difference of this signal to the carrier without modulation would be $\pm5$ radians.
To compute the change in frequency you need to find the compute $\frac{d}{dt}$ of the argument of the cos (divided by $2\pi$) and subtract the carrier frequency.
$$\frac{d}{dt} 10^8\pi\cdot t + 5 sin(2\pi\cdot 10^3 t) = 10^8\pi + 2\cdot 5\cdot 10^3\pi cos(2\pi\cdot 10^3 t)  $$
Divide by $2\pi$ to get Hz instead of radians.
$$Frequency = 5\cdot 10^7  + 5\cdot 10^3 cos(2\pi \cdot 10^3t)$$
Subtract off the carrier of 50MHz and you have a maximum frequency excursion of $\pm 5$ k Hz


How to tell if a high-res flac file has been upsampled from a CD-quality file? Because I might encode my music in lossy format for other devices or edit them, I prefer to keep a high-res copy of my music. However, I know that it is possible to upsample a file using ffmpeg and it would appear the same as an authentic high-res file with ffprobe. Is there anyway I can tell apart a "fake" high-res file (converted from CD-quality or mp3/aac) apart from a "real" high-res file (exported high-res by the producer)?

Let's get some basics out of the way first.
An audio format is basically defined by three different things: sample rate, bit depth and compression type. The compression type can be either lossless (WAV, PCM, FLAC, ALAC) or "lossy" (AAC, MP3, OGG, MQA, etc). Lossy compression works by shaping the quantization noise so it's masked by the remaining content and inaudible. Lossy formats are characterized by their bit-rate: the higher the bit-rate, the less quantization is done.
The typical life cycle is the following: recording, mixing and happens at 96kHz floating point. When that's done, there are typically two masters made: a studio master at 96kHz/24-bit and a CD master at 44.1kHz/16-bit.
The label will send one of these two to the distribution channel, e.g. Spotify. Which one depends on the contract. Spotify will then create the formats that they need for their streaming platform, which almost always includes lossy compression. For example Spotify uses OGG at 320 kb/s. Apple Music uses AAC at 256 kb/s. Which format gets to the customer depends on the contract between the customer and the distribution.

Is there anyway I can tell apart a "fake" high-res file (converted from CD-quality or mp3/aac) apart from a "real" high-res file

Here are a few things to look for in a spectral analysis

Each format has it's own rolloff: CD drops like a rock at 20 kHz and MP3 drops steeply at 16 kHz. If you have a 96 kHz file with these sharp drop offs, it's likely been up-sampled.
Inspect the content above 20 kHz. If there is random "noise like" features in there, it's probably genuine. If it has very little content and/or the content looks like a low-pass filtered mirror image of the content below 20 kHz, it's been up-sampled.
You can look at correlation at high frequencies. For a genuine recording this will mostly be uncorrelated. If there is significant correlation, it's a potential sign of "joint stereo coding" which could hint at lossy compression.
Look at the recording date: If it's been made before 1990 it's almost guaranteed to be up-sampled. There never was a digital studio master and the best they can do is to sample a tape master.

Some formats like MQA, include authentication, i.e. they put a "seal" on the content and can detect whether the "seal" has been tampered with. It's also possible to put a watermark on the content.
Update
You may be able to detect bit padding as well: Any "decent" CD quality file have been properly noise shaped and dithered. It's possible that signs of this will survive the up sampling process can be detected. I've neve tried that, so I don't know how feasible it is, but it sure would be interesting to try.


Reconstruct images from PCA reduced dimensions with NN I was reading this Medium post and I had the idea to reconstruct the original images with a convolutional neural network instead of applying the inverse transform method. The problem is that I don't know how to define the input of the CNN. From what I know, PCA reconstructs the original images as a linear combination of the eigenvectors, so I guess the features should be the k eigenvectors for each original image, but having the same features for different labels doesn't make sense. Is it possible to recontruct the images from the reduced dimensions? If possible, what approach is the one you think would work?

I'm not sure I understand your question. You want to create a neural network able to reconstruct an image thanks to its PCA ? This could be some sort of autoencoder that uses PCA for its latent representation. The encoder of the AE would simply be the PCA (AE isn't actually relevant here, nothing trainable in the "encoder", but I like the parallel since we're talking about generative machine learning/up-scaling 2D representations).

  If possible, what approach is the one you think would work?

Since an image will have 2D eigenvectors of already good shape, you could simply concatenate them, flatten, use a mean, use a concatenation of different 2D statistics/transformations, combine the results of several sub-networks starting from different 2D statistics of the eigenvectors (or the eigenvectors themselves), etc. as an intermediate representation. This is the tricky choice, following which you'll have to deal with typical AE decoder architecture design.

An illustration of the sub-networks idea (replace 1D inputs and 1D output with 2D data in your case):

  Is it possible to recontruct the images from the reduced dimensions?

If you're starting from reduced dimensions, it's not a problem: you can use either the so-called "deconvolution" (transposed convolution is the actual name) for direct 2D upscaling (see here for Keras implementation, and here for the equivalent PyTorch function), or simply if working with a flattened input/intermediate representation, re-shape a 1D large enough neural network layer, to go back to the original image shape.

Finally, note that if you use a neural network as it is commonly defined, your transformations will involve non-linearities.

That being said, not sure it's going to give you a great result, worth a try ! IMO you should try feeding each (or the n-largest) eigenvector into their own very simple sub-network, before combining their respective intermediate representations in a common final network that will output the right shape. Use simple architectures since you probably want to compare with linear combinations of eigenvectors.


Help with deconvolving photon reads distributed across neighboring pixels I'm building a radiation detector that collects photons in CCD pixels and we can relate the energy of the photon from the intensity of the pixel. To test the detector, I took the following spectrum of Iron 55, which is expected to have a single peak at 6 keV:

The tail at lower values of the energy is interesting because it's not supposed to be there. Taking a closer look at pixels in this range, I see that they are clustered next to each other, with what seems to be energies from a single photon distributed into multiple neighboring pixels.

I am intrigued by the possibility of getting a cleaner spectrum by collecting these pixels together. I could try manually going through and clustering at each pixel, but I thought a more elegant way would be to deconvolve the image using a simple point spread function. In the result of this deconvolution I need 2 things to be true:

the energy of the system should be conserved, so the sum of the pixels (the $L_1$ norm) should be conserved.
there should be no negative pixel intensities, since those are unphysical.

I did that, using a Gaussian PSF with a half-pixel standard deviation. However, I find that I either need to violate condition 1 of my requirements, or condition 2 of my requirements. Just using the gaussian PSF with $L_1$ norm 1, which conserves the energy, gives me a plot like this:

The white pixels are negative, unphysical values. At this point I'm stumped, but also obsessed with this problem, because it seems like it should have a clean solution. Any advice?

There is some work on resolving point sources from their smoothed signal called super-resolution, along with some elegant methods to solve it using convex-programming. The idea is that if you know that your signal consists of spikes (single points of energy), you can resolve them from the smoothed signal, even when the smoothing causes interference between spikes.

I used this method in my PhD to find the amplitude of 1D spikes from low-pass filtered signals, using the research published here (not mine). I didn't look at 2D signals like yours, but I recall that some work has been done on it, so you might be able to find some code.

For your method, I'm guessing the negative pixels occur because the smearing isn't exactly modelled by your PSF

What about something simple:

Assume that the smear occurs within a 3 x 3 radius.
Convolve the image with a 3 x 3 filter of all ones. This will add all the values in that area.
Find the points of local maxima in a 3 x 3 radius in the original image.
Replace the value of those points with those of the convolved image and zero the rest.


How is binary data related to double precision data? I am using wavread in MATLAB, which outputs the data in a .wav file in double precision. I am wondering how this relates to binary, such as 0's and 1's. Thank you!

I don't know much about Matlab's wavread, but double precision is a standard for describing a "floating point" value. That is, a value that always has the same number of (binary) digits of precision.
A double precision floating point number in binary looks like this:
Each of those boxes represents a bit (short for "binary digit", or a 1 or a 0). It works like scientific notation: the exponent and fraction are both base-2 representations of numbers, and the number this represents is (with some exceptions for things like 0 and NaN) $$(-1)^{\text{sign}}\times 2^{\text{exponent - exponent bias}} \times 1.\text{fraction}$$
Where the "exponent bias" is there to allow you to represent negative exponents, and the "." in $1.\text{fraction}$ is a "binary point" instead of a decimal point.
You can read about the particulars of the double precision format here, which is where I got the image and formula.


Time-Frquency Resolution issues While Im studying wavelet transform, I have this questions in my mind which can't find it's answer:

I understand that when we take longer time windows and take the fourier transform we would suffer in the time resolution, But Why if we take a shorter time window we would suffer in frequency resolution?
why longer time intervals are   needed for high frequency signals? and shorter time intervals are needed for low frequency signals? In another word, Why we want good frequency resolution for high frequencies and good time resoution for low frequencies?
If I have a signal contains frequencies from 0 Hz to 50 Hz, How can I define the range of low and high frequencies?

I understand that when we take longer time windows and take the fourier transform we would suffer in the time resolution, But Why if we take a shorter time window we would suffer in frequency resolution?

Imagine that your signal is two different sinusoids that are close in frequency and they both start at phase 0.  If the time window is short their end phase will be almost identical because their frequencies are close.  Thus, they are difficult to distinguish with a short time window.

If the time window is long the end phase will be quite different- eventually they will be opposite phases (180 degrees out of phase), and thus very easy to distinguish.  That is why longer time windows give better frequency resolution.

  why longer time intervals are needed for high frequency signals? and shorter time intervals are needed for low frequency signals? In another word, Why we want good frequency resolution for high frequencies and good time resoution for low frequencies?

You do not need longer time intervals for high frequency signals.  You'll need to clarify what you're trying to get at on #2 if you want a better response.

  If I have a signal contains frequencies from 0 Hz to 50 Hz, How can I define the range of low and high frequencies?

Again, I'm not sure what you mean.  What you define as a "high frequency" or "low frequency" is completely arbitrary.  What you will need to do is sample the signal at more than 100 Hz to make sure that the Nyquist frequency is above 50 Hz.


Spatial Aliasing - Wrap Around F-K Spectra I've been using F-K Filter, for a while, but I guess I never had good basic understanding about it Math. Someone asked me what is the cause of frequency wrap around in F-K Spectra plot ? I know it's because of aliasing. But if somebody please elaborate more on the cause of this wrap around? Simple Math explanation perhaps.

Thanks

:)

Aliasing and frequency wrap is a consequence of violation of Nyquist-Shannon sampling theorem which states that a continuous signal must be discretely sampled at least twice the frequency of the highest frequency in the signal.Hence we need to briefly go into the mathematics of the same.

Let x(t) be a continuous signal, y(n) is discrete signal where y(n)= x(nT). Please note that x(t) will have a continuous time Fourier transform (CTFT) while y(n) will have a Discrete Time Fourier Transform (DTFT)

We can construct a mathematical model of sampling using the continuous Dirac delta function p(t) :
$ p(t)= \sum_{k=-\infty}^{k= \infty}\delta(t-KT) \forall t \in Real $

Let  w(t) =x(t)p(t).We can show that CTFT of continuous function w(t) is DTFT of y(n). Writing the frequency domain representation of w(t)

$W(\omega) = 1/2\pi \  X(\omega) * P(\omega) = \int_{-\infty}^{\infty} X(\Omega) P(\Omega -\omega)d\Omega \ \ \  $ , where * is convolution

now CTFT of p(t)  can be written as

$P(\omega) = 2\pi /T \sum_{-\infty}^{\infty} \delta(\omega -k\ 2\pi /T) $

thus $W(\omega) = 1/2\pi \int_{-\infty}^{\infty} X(\Omega)2\pi /T \sum_{-\infty}^{\infty} \delta(\omega -k\ 2\pi /T) d\Omega \\
=   1/T \sum_{-\infty}^{\infty} \int_{-\infty}^{\infty} X(\Omega)  \delta(\omega -k\ 2\pi /T) d\Omega \\
=1/T \sum_{-\infty}^{\infty}X(\omega -k\ 2\pi /T) \   \ \ \ \ \ \ \ \     ...using\  \ shifting\  property  $ ,

Thus we can say that $Y(\omega) = 1/T \sum_{-\infty}^{\infty} X((\omega -2\pi k)/T) $

So we can say that DTFT of y(n): $Y(\omega) $, is a shifted and repeated version of CTFT x(t) i.e. $X(\omega)$. The DTFT is the sum of the CTFT and its copies shifted by multiples of 2π/T.This is shown in the figure bleow. The frequency axis is normalized to −π/T < ω < π/T.

If X(ω) = 0 outside the range −π/T < ω < π/T, i.e. x(t) has no frequency greater than nyquist, then the copies will not overlap the range −π < ω < π and there is no problem as seen in figure above.

However if X has non-zero frequency components higher than π/T (fs/2 or nyquist). There will be overlap causing a wrap around in f-k domain. Refer below:

Notice that in the sampled signal, the frequencies in the vicinity of π are distorted by the overlapping of frequency components above and below π/T in the original signal causing wrap in F-k domain. Now in often spatial sampling is low as in case of Seismic or MRI due to constrains on number of detectors hence there is a aliasing in wavenumber. This spatial aliasing appears as wrap around in Frequency wavenumber Spectra.

Source of Images :R. G. Lyons: Understanding Digital Signal Processing (2nd Edition)

Reference : EECS, University of California Berkeley


Impulse response of non-causal system Suppose I have the following equation, where $x[n]$ is the input and $y[n]$ is the output: $$x[n]-x[n-1]=-10y[n-1]+2y[n]+2y[n+1].$$

The question asks me to compute the impulse response from this equation.

Attempts: $x[n]=1$ when $n=0$; $0$ otherwise. I need to find $h[n]=y[n]$, but, from the equation, it always depends on the future value $y[n+1]$, which is unknown.

One way to solve such a difference equation is to use the $\mathcal{Z}$-transform. Transforming the given difference equation yields

$$X(z)\left(1-z^{-1}\right)=Y(z)\left(-10z^{-1}+2+2z\right)\tag{1}$$

which results in the following transfer function:

$$H(z)=\frac{Y(z)}{X(z)}=\frac{1-z^{-1}}{2z+2-10z^{-1}}=\frac12\frac{z^{-1}-z^{-2}}{1+z^{-1}-5z^{-2}}$$

The poles of $H(z)$, i.e. the zeros of the denominator polynomial in $(1)$, are both real-valued:

$$1+z^{-1}-5z^{-2}=(1-az^{-1})(1-bz^{-1})\tag{2}$$

with $a=\frac{-1+\sqrt{21}}{2}\approx 1.79$ and $b=\frac{-1-\sqrt{21}}{2}\approx -2.79$. Note that $1<|a|<|b|$ is satisfied, so both real-valued poles are outside the unit circle of the complex plane.

With $(2)$, $H(z)$ can be rewritten as

$$H(z)=\frac12\frac{z^{-1}-z^{-2}}{(1-az^{-1})(1-bz^{-1})}=\frac{1}{2(b-a)}\left[\frac{(1-a)z^{-1}}{1-az^{-1}}-\frac{(1-b)z^{-1}}{1-bz^{-1}}\right]\tag{3}$$

This transfer function can have three possible regions of convergence (ROCs), each of them corresponding to a different impulse response:

$$\begin{align}(1)\quad &|z|>|b|>|a|\\
(2)\quad & |b|>|z|>|a|\\
(3)\quad & |b|>|a|>|z|\end{align}$$

ROC $(1)$ corresponds to a right-sided impulse response, ROC $(2)$ corresponds to a two-sided impulse response, and ROC $(3)$ corresponds to a left-sided sequence. Only ROC $(3)$ includes the unit circle $|z|=1$, so only the impulse response corresponding to this ROC results in a stable, yet anti-causal, system. The three impulse responses corresponding to the three ROCs given above are

$$\begin{align}(1)\quad & h_1[n]=\frac{1}{2(b-a)}\left((1-a)a^{n-1}-(1-b)b^{n-1}\right)u[n-1]\\
(2)\quad & h_2[n]=\frac{1}{2(b-a)}\left((1-a)a^{n-1}u[n-1]+(1-b)b^{n-1}u[-n]\right)\\
(3)\quad & h_3[n]=\frac{1}{2(b-a)}\left(-(1-a)a^{n-1}+(1-b)b^{n-1}\right)u[-n]\end{align}$$

These impulse responses describe three different systems, all of them characterized by the same transfer function given by Eq. $(3)$. $h_1[n]$ corresponds to a causal unstable system, $h_2[n]$ to a non-causal unstable system, and $h_3[n]$ to an anti-causal stable system.


What is the impulse response of the system with this difference equation: The difference equation is written:

$$y(k)-\frac{5}{12}y(k-1)+\frac{1}{24}y(k-2)=x(k)-\frac{1}{2}x(k-1) $$

The system is linear, invariant and causal.
For $x(k)=δ(k)$, the equation will be:
$$g(k)-\frac{5}{12}g(k-1)+\frac{1}{24}g(k-2)=δ(k)-\frac{1}{2}δ(k-1)$$

I calculated $g(0)=1$, $g(1)= -\frac{1}{12}$, $g(2)= -\frac{11}{144}$ and $g(3)= -\frac{49}{1728}$ but i can't find any relation between them.

Starting from the difference equation
$$y(k)−\frac{5}{12}y(k−1)+\frac{1}{24}y(k−2)=x(k)−\frac{1}{2}x(k−1),$$
you can express the $Z$-transform on both sides ($X(Z)$ and $Y(Z)$ denote the $Z$-transforms of $x$ and $y$) and get:
$$Y(Z)\left(1−\frac{5}{12}Z^{-1}+\frac{1}{24}Z^{-2}\right)=X(Z)\left(1−\frac{1}{2}Z^{-1}\right).$$

Details on this procedure are quite well explained in $Z$ Transform of Difference Equations.

The transfert function $H(Z)$ is given by the ratio of the output by the input:
$$H(Z) = Y(Z)/X(Z) = \frac{1−\frac{1}{2}Z^{-1}}{1−\frac{5}{12}Z^{-1}+\frac{1}{24}Z^{-2}}.$$
As the denominator is (now) a simple polynomial of degree 2 in $Z^{-1}$, you can quite easily expand the rational fraction $H(Z)$ in $Z^{-1}$ to obtain a series $\sum_{k=0}^{+\infty} h(k)Z^{-k}$ by using formal expressions like:
$$\frac{1}{1-X} = 1+X^{-1}+X^{-2}\ldots$$
and
$$\frac{1}{1+aX+bX^2} = c/(X-r_0)+d/(X-r_1)$$
with $r_0$ and $r_1$ the roots of the degree-2 polynomial $1+aX+bX^2$.

With some computations with the numerator, you should be able to retrieve the solution to your expression, with simple formulas for your $h(k)$, providing you with the impulse response you have looking for.

You can find additional material in Partial fraction decomposition. If you prefer a more self-contained, note-like document, I would suggest the following Difference equation solutions from impulse response, which provides an alternative explanation.


Why do we need FFT pairs for OFDM? this is a basic question, but i barely find any tutorial to explain this.  Most of tutorial dont talk about why use IFFT and FFT, they just say "oh, at transmitter, we just take this information time domain as frequency domain, then we need IFFT". I have no idea why we want to treat as frequency domain when it was time domain. There must be some good reasons, right?

So, why we need the FFT pair for Tx & Rx for OFDM? Why not something else ? Someone could help me with some pictures?  Thank you again, friend.

The short answer is: Complex exponentials are eigenfunctions of an LTI system and the IFFT modulates different complex exponentials with the QAM data.

More elaborate: The wireless channel can (if it is time-invariant for the duration of one OFDM symbol) be modeled as an LTI system with impulse response $h[n]$ (in discrete-time). So, the received signal is the convolution of the input with $h[n]$. At the receiver, you want to know the data, so if you did directly transmit the data, you would need to perform a deconvolution at the receiver.

One can understand the FFT/IFFT to exploit the convolution theorem (convolution in time is multiplication in frequency). We transmit a signal in time domain at the transmitter, and transform it into the frequency domain at the receiver. Then we know, that the frequency response of the channel has multiplied the frequecy domain of the transmit signal. So, it becomes obvious that we directly define the frequency domain of the transmit signal with our QAM symbols.

This is not the whole story (circular convolution with the impulse response must be achieved such that the convolution theorem holds in discrete, and we use a Cyclic prefix to obtain it).

I've written some articles about OFDM, namely an OFDM walkthrough example and about the CP . Maybe these articles can help you understand better. I've always wanted to write something about the math behind OFDM, but didn't find the time yet.


Separating a signal into its constituent components First off, I apologize for possibly incorrect use of jargon, as I am new in this area.

I have a type of signal that I am assuming can be described by an arbitrary number (let's say $n$) of basis vectors, which are linearly independent. I am measuring a signal through one sensor, which gives me a vector of amplitudes. I am able to measure an arbitrary number (let's say $m$) of mixtures of the basis signals, which gives me $m$ measured signals.

The length of my vectors is 150000, sampled at a frequency of 2.5 MHz.

So here is what I think my system of equations will look like.

$$ Y = A X $$

where $Y$ is the matrix containing the mixtures (measurements) and has dimensions of $ m \times 150000 $, $X$ is the matrix containing the basis signals, and has dimensions of $ n \times 150000 $, and $A$ is a coefficient matrix that has dimensions of $ m \times n $.

I need to solve for $A$. However, since $X$ is non-square ($ n \ll 150000$), it does not have an inverse. I guess I could calculate the pseudo-inverse. Does this make sense?

I have also looked into Principal Component Analysis and how it is used in various types of spectroscopy. It seems to me that PCA assumes a set of basis signals, and I don't know whether that has a physical meaning in my case, since my signals are stochastic (but stationary). I used the PCA routine on matrix $Y$ in MATLAB, and it gave me some signals very similar to the input ones, but I did not know what to do with them.

Any ideas?

I think I do have a valid physical model for my signals. Should I perhaps use model-based blind source separation (BSS) ?

Thanks in advance for taking the time.

You do a QR decomposition of $X^T$. Then the system reads as $$Y=AR^TQ^T.$$ Now multiply with $Q$ from the right, $$YQ=AR^T,$$ which already should give you hints on the reliability of the eventual solution, since ideally, the block of $YQ$ right of the leading $m\times n$ block should be zero. If that is true, then $$A=[YQ]_{[1:m,1:n]}[R^T]_{[1:n,1:n]}^{-1}$$.

Note that in the QR decomposition using Householder reflections, you can apply the reflections computed for $X$ simultaneously to $X$ and $Y$, if $$S=I-2 vv^T,$$ then $$X_{next}=XS=X-2(Xv)v^T\text{ and }Y_{next}=YS=Y-2(Yv)v^T$$ and $$X_{final}=R^T,\; Y_{final}=YQ$$.


Power spectral density: Why are these two methods equal? The power spectral density can be calculated in two ways:

by doing the Fourier_transform of the autocorrelation
by doing (abs(X(f)).^2 where X(f)=fft(x(t))

Can you explain me passage by passage why these are equal?

I'll try to explain better my problem, because I tried to simplify it, but maybe (since I am not an expert) I have worsen the situation.

I had to do the power spectral density of a random signal (the velocity of a CFD simulation).

Firstly I did the FFT of the autocorrelation:

   v is my signal
   r =xcorr(v,v)
   PSD = abs(fftshift(fft(v)))

After that, I realised that in order to be verified the fact that the integral of the PSD is equal to the rms^2 of the signal, my result has to be scaled by the factor 1/(fs*N), as it is shown in PSD, Autocorrelation, and Noise in MATLAB.

Unfortunately I did not understand why this factor was necessary and I did not find any explanation about this online or in my books. So I looked at this video that gives an explanation using the Parseval Theorem, How to scale FFT and PSD in MATLAB.

So I thought if I demonstrate that abs(fftshift(fft(v))) is equal to abs(fft(v)).^2 then I can demonstrate why I need that scaling factor 1/(fs*N) in order to satisfy my check.

I hope you can forgive me if I did not explained everything from the beginning, but I tried to make it simple (I obviously failed). I hope you can help me. I am struggling with this demonstration for one week and none can help me.

Regarding what we discussed in the comments, here is an explanation of the different concepts I was referring to:

The concrete usage of terms may vary, but my interpretation is this: For a stochastic process $x(t)$, the autocorrelation function is defined as $$r(\tau) = \mathbb{E}\{x(t) \cdot x(t+\tau)\}.$$ Note that the process needs to be stationary, otherwise the autocorrelation function depends on $t$ as well. Its power spectral density is given by $$R(f) = \mathbb{E}\{|X(f)|^2\}.$$ These two quantities are linked to each other via the Wiener-Khinchin-Theorem which shows that $$R(f) = \int_{-\infty}^{\infty} r(\tau) {\rm e}^{-\jmath 2\pi f \tau} {\rm d} \tau,$$ i.e., they are Fourier transform pairs.

Since the expectation is hard to implement in practice (we typically do not have an infinite ensemble available), we need to resort to other methods. When the process is ergodic, ensemble average can be replaced by temporal average and we may, e.g., estimate the autocorrelation function from a window $T$ via $$r_T(\tau) = \frac{1}{T} \int_T x(t) \cdot x(t+\tau){\rm d}t.$$ Note that $r(\tau) = \lim_{T\rightarrow \infty} r_T(\tau)$. On the other hand, considering using $R_T(f) = |X_T(f)|^2$ as an estimate for the power spectral density (referred to as the periodigram), where $X_T(f)$ is the spectum of a $T$-windowed signal leads to an estimate that does not converge since $S(f) \neq \lim_{T\rightarrow \infty} S_T(f)$.

Now, I'm not sure what exactly is the question. Is it the explanation of the Wiener-Khinchin-Theorem for the stochastic signals? Or is it the deterministic equivalent as, e.g., discussed in this answer?

*edit: Regarding your clarification that you posted as a reply to your question(?): My feeling is this question has little to do with ACF and PSD being Fourier transform pairs, instead you seem to be worried how to normalize the result of a Discrete Fourier transform (carried out via an FFT) such that its units make sense physically. It's probably worth asking this in a separate question, but the short version is here. There are two things to be considered:

There are different conventions how to normalize DFT and IDFT. Let $\mathbf{F}$ be a DFT matrix with elements ${\rm e}^{-\jmath \frac{mn}{N}}$, then you can either (i) define DFT via $\mathbf{D} = \mathbf{F} \cdot \mathbf{d}$ in which case you need to define the IDFT via $\mathbf{d} = {\color{red}{ \frac{1}{N} }} \mathbf{F}^H \cdot \mathbf{D}$ or (ii) define the DFT via $\mathbf{D} = {\color{red}{ \frac{1}{N} }}\mathbf{F} \cdot \mathbf{d}$ which leads to the IDFT $\mathbf{d} =  \mathbf{F}^H \cdot \mathbf{D}$ or alternatively (iii) you let the DFT be $\mathbf{D} = {\color{red}{ \frac{1}{\sqrt{N}} }}\mathbf{F} \cdot \mathbf{d}$ which leads to the IDFT $\mathbf{D} = {\color{red}{ \frac{1}{\sqrt{N}} }} \mathbf{F}^H \cdot \mathbf{d}$. Different conventions are in place, Matlab uses the first by default, whereas Maple will prefer the third.
The Fourier transform changes the units: the ACF has a unit of power (say, Watt), the Fourier transform is an integral over time, hence it adds seconds, so that the unit is Watt seconds or better Watt/Hz. It's a density anyways (how do Watts distribute over Hertz?). If you want to approximate a Fourier transform via a Discrete Fourier transform, you have to be super careful, since the DFT/IDFT does not change the units. If your input to the DFT is a sequence with unit, say, Watt, its DFT coefficients will have the same unit. If you want to mirror what happens in a Fourier transform, you need to multiply the result with the sampling period of the samples in time domain. In other words: If you have a function function $u(t)$ with spectrum $U(f)$ and you want to approximate $U(f)$ with a DFT, you need to consider a finite window (say, length $T$) and a finite sampling interval (say, $t_0 = T/N$). Then, ${\color{red}{t_0}}\cdot D[\mu]$ is an approximation for $U(\mu f_0)$ where $f_0 = 1/T$. It's an approximation since sampling as well as truncation can cause errors (one of them will for sure).

I think the factor you were looking for can very likely be a combination of these two effects, your $1/(f_s N)$ is the same as my $t_0/N$.


How to identify causality, stability and ROC from the pole-zero plot? To preface, this is not a homework related question but purely for self-study purposes.

If I am given the following Pole-Zero Plot: (Source: Berkeley Exam1)

How would I go about trying to determine whether the system is causal or stable?
Similarly, why should the region of convergence be $0.5 < |z| < 2 $ for stability?

I know that for stability, the ROC must include the unit circle.  And for causality, $h[n]=0$ for $ n<0$.  Other than that, I am having trouble applying those rules to a pole-zero plot.

For systems with more than one pole (with different radii), there are $3$ types of ROCs:

inside the circle with a radius corresponding to the smallest pole radius; the corresponding time-domain sequence is left-sided
outside the circle with a radius corresponding to the largest pole radius; the corresponding time-domain sequence is right-sided
in a ring defined by the radii $r_1$ and $r_2$ of any two poles with $r_1<r_2$, with no other pole at a radius in the interval $(r_1,r_2)$; the corresponding time-domain sequence is two-sided

In your example you have case $3$ above, because that's the only ROC that includes the unit circle, and, consequently, corresponds to a stable system. Both other possible ROCs (inside the circle with radius $0.5$ and outside the circle with radius $2$) do not include the unit circle, and, consequently, do not correspond to impulse responses of stable systems. Since a ROC equal to a ring between two poles corresponds to a two-sided sequence, the system can't be causal. A causal system has a ROC outside a circle (case $2$ in the list above).


What's the expression of the noise for a RF signal? I am confused about how to give the expression of noise. Can someone tell me if my equations are correct?

If the received RF signal is expressed as:

where

If I want to add a AWGN, can I give the expression of noise as the following?

If not, what is the correct equation of the noise? Also, in the MATLAB simulation, I often add noise to a baseband signal by using: sqrt(NoisePower) * randn(). What should I do in the case of bandpass signal? Thanks guys!

Yes, you can represent your bandpass noise as:

$ n(t) =n_I(t)\cos(2\pi f_c t)-n_Q(t)\sin(2\pi f_c t),$

where $n_I(t)$ and $n_Q(t)$ represent the inphase and quadrature components of the noise signal. In matlab you can generate this sequence as

n = sqrt(Npwr/2) *(randn(1,len)+j*randn(1,len))

The scaling factor is different than in the baseband case because both the inphase and quadrature components have unit power (before scaling), so when added together the total power is equal to 2 and therefore the $\sqrt{2}$ factor is needed in the scaling factor to get unit power.


differences in fundamental frequency/pitch between music and speech I am developping a program, that extracts very expressive frames from a video (meaning: frames, which could summarize a video) for a university seminar. One of my ideas was, to analyze the audio data (e.g.: if it's very loud, there must be sth. happening and so on). I'm using Praat for this analysis. Now, I stumbled upon the pitch, which is supposed to be the fundamental frequency (formant f0). Can anyone tell me what usually happens to the pitch, when somebody talks or when you can hear only music? And do you have any other ideas what the audio data of a film might be able to tell me?

I don't think pitch information is relevant for what you want to do.

The variation of pitch during speech is known as intonation, and can convey emotions, indicate if a sentence is a question etc. However, there is no universal rule as to how pitch variation patterns are mapped to meaning - this is quite language dependent ; and some languages sound "flatter" in terms of pitch contour than others.

While there are some exceptions (say a solo flute performance), music is polyphonic in that several instruments are playing at the same time - and to make things worse some of those instruments might be playing several notes at the same time. As a result, at a given point in time, music is unlikely to contain just a single fundamental frequency ; and when applied to music, Praat's f0 estimator is very likely to yield a noisy and unreliable pitch contour, sometimes jumping from one instrument to the other.

Regarding automatically annotating movies from the audio, here are a few ideas:

Looking for sharp increases in loudness to detect important events as you suggested.
Use speech/music discrimination to help detect scene boundaries.
Train classifiers for common sound effects (explosions, gunshots, punches), and use this to detect action/violence scenes.
Use speaker diarization techniques (along with some face recognition) to detect all lines spoken by a particular actor.


Why is the first axis of STFT frequency bin When taking the STFT of an audio signal, the first axis is the frequency bins by convention, and the second axis is how it elapses through time. Is there a reason?

It is really a matter of convention of the communities that use STFT. Medical devices tend to have time scrolling left to right or over writing left to right. For passive SONAR, most displays that Im familiar with have time water falling down and frequency increasing from left to right. There are exceptions as well. I've also seen amplitudes as marching mountains, color coded, and intensity coded as well. Much depends on the amount of time history that should be displayed and the time duration of signal events that are of interest.  There are also perceptual consideration based on human vision that come into play and effect how STFT data is displayed. There are also places where an STFT display is the input to image processing algorithm.


why is the DFS of a delta function equal to 1 I have a x[n] = $\delta$[n].
By formula is should be

$$

X[k]= \sum_{n=0}^{N-1} \delta[n]W_N^{kn}

X[k]= \sum_{n=0}^{N-1} e^{-j2*pi*kn/N}

$$

The formulae isn't showing for some reason. I took a screenshot of what I got here: http://imgur.com/6j0Ibgu

basically why is the summation of an exponential term going to 1. $W_N=e^{-j*2*pi*k*n/N}$ in this case. I tried to prove it via $\sum\alpha^k=\frac{1-\alpha^N}{1-\alpha}$ but it doesn't work for me.

$X[k]= \sum_{n=0}^{N-1} \delta[n]W_N^{kn} \quad$, where $\>W_N^{kn}=e^{−j\,2\pi k\,n/N}$

HINT:

What is the value of $\delta[n]$ when $n \neq 0 \>$ ?


Calculating an incoherence property from sub-optimal sampling patterns EDIT (after comments and subject matter review)

CS is based on a choice of a sensing basis $\Phi$ relative to a representation basis $\Psi$. Using an "Incoherence Property" $\mu$ that measures the correlation between vectors from each basis it has been shown that the CS paradigm can achieve good reconstructions if the incoherence between $\Phi$ and $\Psi$ is low. Examples of low $\mu$ come from: The spike basis ($\Phi$) and the Fourier basis ($\Psi$), the noiselet basis ($\Phi$) and the wavelet basis ($\Psi$), a random matrix ($\Phi$) and a fixed basis ($\Psi$).

updated question is:

Is there a way to calculate the $\mu$ when the user is not able to chose $\Phi$ ahead of time to sample the signal, i.e. when the samples are already provided through some arbitrary sampling scheme that is sub-optimal vis-vis CS. My aim is to show/quantify, via $\mu$, that the arbitrary sampling scheme I have been given is sub-optimal if one is expecting CS to work.

ORIGINAL

With respect to Matrix Completion and Compressive Sampling (CS) I'm trying to understand how to calculate an incoherence property μ between two bases Φ and Ψ. Getting this incoherence is important because if Φ and Ψ are highly correlated there is little chance of succesfully reconstructing a signal from sparse samples. It is stated here (page 3) that μ is given by:

where n represents the number of elements in a matrix M - say an image signal.

I understand Φ to be a sensing basis and Ψ to be a sparse representation basis. I am using uniform random sampling to get a set of a samples from M. But I would also like to try arbitrary sampling patterns. (I'm not using wavelets or noiselets or Fourier coefficients - although I would like to try this eventually)

My question is:

How do I actually obtain Φ and Ψ?

$\Phi$ is the matrix that represents the way you sample your signal $x$. Actually,  $\Phi$ can be an identity matrix with some rows eliminated, which means you are picking a subset of $x$.

$\Psi$ represents the basis that you choose to expand $x$:
$$x = \Psi c$$
A very simple example is that $x$ is a cosine function (thus containing only one frequency component), then you can take $\Psi$ to be the Discrete Fourier (or discrete cosine) matrix and $c$ is just the coefficient of the specific frequency, and it is pretty sparse (all the coefficients on other frequencies will be zero).

Then the compressed signal $b$ can be represented as:
$$b = \Phi \Psi c$$

Actually you can also view $\Phi \Psi$ as selection of the random rows of $\Psi$, then
one possible option of $\Phi$ and $\Psi$ is (in Matlab):

Psi = fft(eye(n));
p = randperm(n);
Phsi = Psi(p(1:m),:); % choose the first m rows from the permutation

Psi is $\Psi$, and Phsi is $\Phi \Psi$.

For your case, yes you can use singular value decomposition on the image space. But I would prefer to

[u,s,v]=svd(M'*M); % spectrogram
Psi = s(1:m,1:m) * v(1:m,:);
U = u(1:m,:);

which indicates that you select the m most magnitude on spectrum components. Columns in v' are the orthogonal vectors that map the frequency component with the corresponding magnitude to spectrum space. Please note that you don't need u as the basis, since it is always equivalent between u * s * v' * x = b, and s * v' * x = u' * b. As a result, the compressed form is converted to:
$$U'b = \Phi \Psi c$$

Hope it helps. Thanks


Triangle-kernel filtering problems at non-integer ratios I'm trying to implement image downscaling using a triangle filter (it's narrow therefore fast, and goes hand in hand with linear interpolation when upscaling), and the triangle function used for filtering has a span proportional to the downscaling ratio.

Alas when the ratio is non-integer, grid artifacts appear (as shown above) on the downscaled image. The following plot shows what happens with the way triangles for each pixel of the input image overlap and add up in those cases:

I don't know what I should do to fix/work around this, but I'm sure this problem must have been encountered by others before (Google was no help).

I went with Jason R's idea of adjusting the linear filtering, except I didn't like the idea of doing a division on every output pixel, using more than just a simple kernel or how normalising that way would ruin the smooth borders. So instead I modified the kernel itself, which it turns out is done by flattening the peak of the triangle.

It seems like a pretty good solution, and you can see how the sum of the triangles progresses into something resembling a Gaussian function.

Here's how it's done, first of all there's a 'knee' joining the flat top with the rest of the triangle. It ping-pongs on the horizontal axis between 0 and 0.5 like so. Note: $n$ is the scale multiplier for the width of the kernel and the size of the downscaling, it's always 1 or lower, it's the reciprocal of the number displayed at the top of the animation. The smaller $n$, the smaller the output image, the wider the kernel over the original image.

$knee = 0.5 - |(n^{-1} \mod 1) - 0.5|$

Next we need to calculate the height of the top of the flattened triangle, which is done like this:

$m = \lfloor n^{-1} + 0.5\rfloor$

$top = \frac{1 - knee \cdot n}{2m - nm^2}$

And finally we turn it into an affine function which we can use directly when calculating weights:

$c_1 = \frac{-top}{n^{-1}-knee}$

$c_0 = \frac{-c_1}{n}$

The above is calculated only once and $knee$, $top$, $c_1$ and $c_0$ are stored and used when computing weights like this:

$$f(x) =
  \begin{cases}
      \hfill top   \hfill & \text{ if $|x| \leq knee$ } \\
      \hfill c_1 x + c_0 \hfill & \text{ if $knee < |x| \leq n^{-1}$ } \\
      \hfill 0   \hfill & \text{ otherwise } \\
  \end{cases}$$


When deriving the power spectral density of stochastic processes, why does taking an expectation allow the $T\rightarrow\infty$ limit to be taken? I am following the arguments presented in the paper AN-255 Power Spectra Estimation, from Texas Instruments, to learn how to derive the power spectral density for a stationary stochastic process, and do not understand one of their steps. I am looking at Page 4 only of the the document for the purposes of this question.

The derivation is done by considering only a single realisation (sample function), denoted by $x(t)$, of the stochastic process. Then they do the usual thing of only considering a truncated version of this signal $x_T(t)$, which allows the Fourier transform of it to exist, denoted by $X_T(f)$ (because one of the requirements for FT to exist is absolute integrability, which isn't satisfied for random signals).

I am confused with how they justify taking the limit $T\rightarrow\infty$ in the right hand column of page 4, shown here:

$\hskip2in$

$\hskip2in$

I understand that you cannot take the limit $T\rightarrow\infty$ directly on the truncated sample function's Fourier transform $X_T(f)$, because it is not defined in that limit (indeed, that was the reason for doing the truncation in the first place).

I also understand that $X_T(f)$ results from a single realisation only of the stochastic process, and if you were to repeat the experiment again you would obtain a different $X_T(f)$. As such, this function is actually a random variable itself (at each frequency $f$), and therefore it makes sense to do an ensemble average over many realisations, and consequently motivates taking an expectation.

What I don't understand is why the expectation solves the problem, and allows you to take the limit $T\rightarrow\infty$ after you have taken the expectation:

$\hskip2in$

What is it about the expectation of $X_T(f)$ that means the limit $T\rightarrow\infty$ can then be taken? I'm having a hard time seeing it... Could someone spell this out plainly for me?

It is not the expectation operator that makes sure that the limit exists. The expectation just results in an ensemble average, which we need to obtain a deterministic function $S(f)$ for the power spectrum.

Assume we're given a deterministic power signal $x(t)$, i.e., a signal with finite non-zero power, and, consequently, infinite energy. Its Fourier transform generally doesn't exist. We can define a truncated Fourier transform

$$X_T(\omega)=\int_{-T}^Tx(t)e^{-j\omega t}dt\tag{1}$$

By assumption, $\lim_{T\to\infty}X_T(\omega)$ doesn't exist. However, the limit

$$\lim_{T\to\infty}\frac{\left|X_T(\omega)\right|^2}{2T}\tag{2}$$

exists because of the finite power assumption, and this is also the way the power spectrum is defined for such deterministic power signals.

If $x(t)$ is modeled as a random signal, we only need to modify $(2)$ by taking the expectation of the numerator to obtain the power spectrum of $x(t)$.


Merits of "Zero-Phase" Zero Padding I understand the reasons for zero-padding a finite signal before FFT.

However, I was reading this book (http://www.dsprelated.com/dspbooks/sasp/Zero_Phase_Zero_Padding.html) where the topic is Zero Phase Zero Padding.

Top image: windowed finite signal.
Bottom image: Zero Phase Zero Padded version of the same signal.

The merit of Zero Phase Zero Padding over "regular" zero padding (whether adding zeros before, after or at both ends of the finite signal) isn't obvious to me. What is its advantage?

EDIT:

To help anyone else who might have this question in the future, I am adding my comments here and things I tried in order to understand this more - hopefully this will save time for anyone with this question in the future.

After @matt-l's comment, I wrote a matlab script for seeing the difference in the magnitude and phase responses of a symmetric signal that was zero-padded in 3 ways (signal000000, 000signal000 and nal000000sig). The third one is Zero-Phase Zero Padding.

The imaginary part of the frequency response in case 3 is very close to 0 (on the order of magnitude of the noise floor - eps in Matlab). However, the unwrapped phase still exists - which, as Matt pointed out might be due to numerical issues.

High-res

The script I wrote is:

NSIG = 128; % Signal length
NFFT = 512; % FFT length

% Make signal
n = 0:NSIG - 1;
w1 = pi/8;
w2 = 3*pi/4; % w1 and w2 are normalized frequencies of the signal

signal = 0.5*cos(w1*n) + cos(w2*n);

% Make signal symmetric about NSIG/2
signal = [fliplr(signal), signal(2:end)]; % Make signal symmetric
NSIG = length(signal);

% Window the input
window = ones (size(signal));
% window = hanning(NSIG)';
raw_in = signal.*window;

%% A. Zero pad at the end only
NZERO = NFFT-NSIG;
in = [raw_in, zeros(1,NFFT-NSIG)];
show_spec1 = fft (in);

%% B. Zero pad on both sides
in = [zeros(1, floor(NZERO/2)), raw_in, zeros(1,ceil(NZERO/2))];
show_spec2 = fft (in);

%% C. Zero phase zero pad
in = [raw_in(floor(NSIG/2)+1:NSIG), zeros(1,NZERO), raw_in(1:floor(NSIG/2))];
show_spec3 = fft (in);

% Prepare to plot
fs = 500;
fres = fs/NFFT;     % Freq resolution
w = fres.*([0:NFFT/2, -NFFT/2+1:-1]);
figure;

%% PLOTS
titleStr = 'A. - s i g n a l 0 0 0 0 0 0 -';
h(1) = subplot (3,2,1); stem (w, abs(show_spec1),'.-'); axis tight; title (titleStr); grid on;
xlabel('Frequency (Hz)'); ylabel('Magnitude');
h(2) = subplot (3,2,2); plot (w, unwrap(angle(show_spec1)),'r'); axis tight; title (titleStr); grid on;
xlabel('Frequency (Hz)'); ylabel('Unwrapped Phase');

titleStr = 'B. - 0 0 0 s i g n a l 0 0 0 -';
h(3) = subplot (3,2,3); stem (w, abs(show_spec2),'.-'); axis tight; title (titleStr); grid on;
xlabel('Frequency (Hz)'); ylabel('Magnitude');
h(4) = subplot (3,2,4); plot (w, unwrap(angle(show_spec2)),'r'); axis tight; title (titleStr); grid on;
xlabel('Frequency (Hz)'); ylabel('Unwrapped Phase');

titleStr = 'C. - n a l 0 0 0 0 0 0 s i g -';
h(5) = subplot (3,2,5); stem (w, abs(show_spec3),'.-'); axis tight; title (titleStr); grid on;
xlabel('Frequency (Hz)'); ylabel('Magnitude');
h(6) = subplot (3,2,6); plot (w, unwrap(angle(show_spec3)),'r'); axis tight; title (titleStr); grid on;
xlabel('Frequency (Hz)'); ylabel('Unwrapped Phase');

This is just about obtaining a symmetric signal after zero-padding. Take a symmetric signal (w.r.t. to time index $n=0$) and append zeros. Due to the implicit periodicity of the time signal used as input to the DFT, it does not matter if you append all zeros to the right, to the left, or on both ends of the symmetric signal. In any case, since the DFT/FFT assumes the first sample in the buffer to be the time origin, the input vector to the FFT needs to be arranged as

[right half of the signal | left half of the signal]

which is basically one period of the periodically continued zero-padded signal between $n=0$ and $n=N-1$, where $N$ is the FFT length, i.e. the length of the signal after zero-padding. Note that now all the zeros end up in the middle of the vector to be transformed. If done correctly, the result of the FFT will be purely real-valued (apart from numerical errors). If the signal is asymmetrical, its FFT will be purely imaginary.


How to obtain mean Power Spectrum from two audio channels? I have a source of noise in a duct, on either side of the noise source I have a microphone, which register the pressure, so that I get two pressure signals: $p_1$ and $p_2$. I want to get the acoustic Sound Pressure Level (SPL) spectrum which can be estimated using a periodogram approach, more specifically using the welch procedure. This function is readily available in Matlab as

pwelch

or in python using

scipy.signal.welch

The Question

The question is, since I have two channels of audio, how do I combine them to get a representation of the total SPL spectrum? Or is the mean SPL more appropriate?

Some ideas:

I know that for two sound sources giving rise to $p_1$ and $p_2$ then we have two cases:

Correlated sound: $<p_{tot}^2> = \frac{1}{2}(p_1^2+p_2^2)$
Uncorrelated sound: $<p_{tot}^2> = \frac{1}{2}(p_1+p_2)^2$

In my case I'm guessing I should consider the time series $p_1$ and $p_2$ to be correlated since they come from the same source but are measured on separate sides, and this would imply that the total SPL spectrum would be given by averaging the resulting spectra given by welch on the individual channels. Thoughts?

The sound pressure level (and derived metrics such as the SPL spectrum) are dependent on the position of the listener, since it is a measure of the pressure of the wavefront moving through space. As such it doesn't make all that much sense to try to combine SPL's from different positions into a single value, except perhaps if the listeners were exactly the same distance from a single omnidirectional audio source and there was no possibility of reflections.

The equations you have posted are for two audio sources, not two listeners, so I am not sure they are relevant.

The only value of any use I can imagine you being able to calculate from two listeners would be the peak SPL / spectrum, this would simply be the maximum value of both channels for any given sample.


What are the basic categories of operations performed on audio signals to create interesting sounds? I don't know if this question will make sense as I'm very new to dsp.

In my limited foray into audio dsp I've come across:

Frequency Modulation
Amplitude Modulation
Additive Synthesis
Subtractive Synthesis

My question is: are these the main fundamental categories of signal manipulation (specifically for audio)?

In other words can all the effects and fancy things that FL Studio with a full suite of plugins can do, be broken down into a series and combination of the above 4 operations? Or are there a whole bunch more? Can a modern DAW with plugins theoretically be broken down into any categories in this manner? Do the above 4 categories even make sense in themselves?!

A bit of background: Essentially I'm trying to make a very basic (but very flexible) audio synthesizer/editor as a project to learn both programming and audio dsp. Using java I started by mimicking basic real world synthersizer-ish classes, calling things like Oscillator and LFO and the like, but every time I learn something new, I find myself having to rewrite everything to fit the new concept into the program.

I'm trying to start again now and come up with the basic objects and structure of the program. In order to that properly I need to know how I should allow these objects to interact etc...

Thanks!

EDIT***

Thanks for the useful comments and answers. I realise I've probably been underestimating the dsp side of stuff and I won't be able to just create a basic template & structure for a simple synthesizer and then just expand it/build on it. I will probably have to repeatedly rebuild again & again as I learn new stuff & the "proper" way to do things... It will take longer, but maybe I'll learn more? Also sorry I don't seem to have enough rep to upvote anyone yet... as soon as I can, I will.

EDIT*EDIT***

I should add that a bit of googling revealed this basic guide to "types of synthesis" which I found relevant and interesting.

The AM / FM category is a bit strange for effects - true, you can get some interesting effects by modulating the amplitude of an input signal ; but what does "frequency modulation" mean for a complex input signal for which you don't even have an accurate frequency representation? You could very well say that every effect is an amplitude modulation, by the $\frac{out(t)}{in(t)}$ signal, but this is not going to help you implement it!

There is no magic single-purpose engine at the heart of all the audio effects provided by music production software. But if you look at the source code of a large suite of audio effects from a DSP developer, here are some "building blocks" that will have been factored out in their own classes and which are shared by many different effects. I don't imply that this forms a good basis to categorize effects ; but from an implementation point of view, those blocks are essential:

Biquad filters.
Fractional delay lines.
Convolution engine, for fast convolution of an input signal with arbitrary-sized impulse responses.
Waveshaper (application of a non-linear function to the input signal in the time-domaine).
Synthesizer building blocks: oscillators, LFOs, ADSR envelopes.
Signal detectors: envelope follower, f0 extractor.

With those blocks you could build:

Synth-like filters or wah-wah: one or two biquads
Auto-wah: envelope follower + envelope + biquad
Flanger / Phaser: LFO + delay line
Chorus: LFO + delay lines
Algorithmic reverb: delay lines (array of parallel and serial comb filters)
Convolution reverb: convolution
Spatializer: convolution (with HRTFs impulse responses)
Guitar amp simulation: convolution, waveshappers
Distortion, Fuzz: gain + waveshapper
Bitcrusher (quantizer): a particular case of waveshapper
Ring modulator with a sine carrier: oscillator
Noise Gate and other dynamics effects: envelope follower
Tape simulation: convolution, waveshapper, envelope follower (for the dynamics processing)
EQ: biquads
Tremolo: LFO
Leslie emulation: convolution + LFO
Vocoder: biquads + envelope follower

While some of these effects are just a matter of patching the different blocks (a flanger truly is a LFO modulating a delay line), some other effects need more complex glue between the parts, that is specific enough not to be factored into a building block of its own.

This is, overall, an interesting set of effects that will cover a lot of ground for music production, but it's not exhaustive and there are indeed some effects which do not fit this framework... Some examples:
* Dirty sample rate reducer (in a bitcrusher): true, it's a multiplication by a dirac comb and then a convolution by a rectangular window... but it's easier to write it as something that duplicate the value of one sample over the N following samples instead of patching an amplitude modulation and a convolution!
* Pitch alteration effects (pitch-shifting, auto-tune) do not fit this framework well. They need a more complex representation of the audio signal to be performed (phase vocoder for frequency-domain methods ; pitch detection and pitch marks for time-domain methods).

I suggest you to play with "modular style" software tools like Pd, Reaktor, Plogue, SyntheEdit... - and try implement effects from the basic building blocks they provide you.