Fixed-point and ARM CMSIS-DSP FIR filters I am using Q1.15 fixed-point values to represent a 12-bit signed integer value. I was wondering if it matters if I use the  11 MSB or LSB of the fractional part to represent my 11 magnitude bits? Only thing I can think of is using LSB might in some cases reduces the chance of overflow errors and maybe allow me to use larger fixed-point FIR coefficients.

But I am using the ARM CMSIS-DSP FIR filtering functions and from what I understand the Q15 set of functions expects Q15 coefficients, Q15 state variables, Q15 result and of course Q15 input data. So in this case I will not be able to increase the precision on my Q15 filter coefficients.

Is there then any difference?

I don't think either will increase your precision as such, although there are benefits to each approach.

Using the lower 11 bits of 1.15 would provide you with 24dB of extra headroom (obviously being careful to sign extend properly). Alternatively, using the upper 11 bits could potentially lower the power in your quantisation error after 1.15 * 1.15 multiplication (for example).

The design is a compromise being headroom and dynamic range, and hence the best route probably depends on the magnitude of your coefficients.

The following example shows how extra headroom can be useful:

Let's say you have an FIR filter (1.15 fractional binary) with 2 coefficients, which are both 0.8 (for some reason). If your input was 0.9 for two samples in a row, your output would be (0.8 * 0.9) + (0.8 * 0.9) = 1.44. This value would lead to an overflow or hard clipping since it doesn't fit into 1.15 (fractional). To avoid this you could first shift the input down by 1 bit, to give yourself 6dB of headroom. This would change the sum to be (0.4 * 0.9) + (0.4 * 0.9) = 0.72, which now safely fits into a 1.15 container.

In a normal system where you're using all 16 bits you'd then have to worry about (or just remember!) the fact that your signal is halved in amplitude. In your case you don't need to worry about this, as your full scale signal is only 12 bits to start with.


Matlab: How to apply constant modulus algorithm in wireless communication I am trying to apply the Constant Modulus Algorithm which is a blind equalization algorithm in CDMA communication. I am following the code for CMA and a great explanation given in link MATLAB : Proper estimation of weights and how to calculate MSE for QPSK signal for Constant Modulus Algorithm

For different signal-to-noise ratio, I want to see the performance of CMA. But I cannot understand where to plug in the CMA part in CDMA. CMA should be at the receiver, but it is not clear to  me exactly how to do so and if I need to apply the same FIR channel that is given in the CMA code. It will be very helpful is somebody can inform / show how to apply CMA. Thank you.
The code for CDMA is

clc;
clear all;
code_length = input('Enter Code Length = ');
data_stream = input('Enter Data Bit Stream = ');
%Generating Walsh Code
code = [-1 -1; -1 +1];
[r1 c1] = size(code);
disp('Walsh code generated = ');
while r1<code_length
    code = [code,code;code,-1*code];
    [r1 c1] = size(code);
    disp(code);
end
code_length = length(code); %length (number of bits) of code
Y = size(data_stream);
N = Y(1); %number of unique senders / bit streams
I = Y(2); %number of bits per stream
T = []; %sum of all transmitted and encoded data on channel
data_received = []; %vector of reconstructed bits at receiver
%show data bits and codes
disp('Data bit stream to be transmitted = '),disp(data_stream);
disp('Walsh codes used for transmission = '),disp(code);
%encode bits and transmit
G = zeros(I,code_length);
for n = 1:N
    Z = zeros(I,code_length);
    for i = 1:I
        for m = 1:code_length
            Z(i,m) = [data_stream(n,i)*code(n,m)];
        end
    end
    G = G + Z;
end
%show channel traffic
for i = 1:I
    T = [ T G(i,:) ];
end
disp('Resulting traffic on the channel = '),disp(T);
%decode and reconstruct
for n = 1:N
    tot = zeros(1,I);
    R   = zeros(I,code_length);
    for i = 1:I
        for m = 1:code_length
            R(i,m) = G(i,m) * code (n,m);
            tot(i) = tot(i) + R (i,m);
         end
    end
    data_received = [data_received ; tot / code_length];
end
disp('Reconstructed data at the receiver = '),disp(data_received);

Output:

Enter Code Length = 4
Enter Data Bit Stream = [ 1 -1  1 -1  1  1 -1 -1 ;
        -1 -1  1  1  1 -1 -1  1 ;
         1  1 -1 -1 -1  1  1 -1 ;
         1  1  1  1 -1 -1 -1 -1 ];

Walsh code generated =
    -1    -1    -1    -1
    -1     1    -1     1
    -1    -1     1     1
    -1     1     1    -1

Data bit stream to be transmitted =
     1    -1     1    -1     1     1    -1    -1
    -1    -1     1     1     1    -1    -1     1
     1     1    -1    -1    -1     1     1    -1
     1     1     1     1    -1    -1    -1    -1

Walsh codes used for transmission =
    -1    -1    -1    -1
    -1     1    -1     1
    -1    -1     1     1
    -1     1     1    -1

MATLAB : Proper estimation of weights and how to calculate MSE for QPSK signal for Constant Modulus Algorithm

Ok, your code is very difficult to understand, but I think that we can answer your question.

There is a group of processing blocks that every communication systems have. Normally, they appear in the following order in almost every simulation:

Bit Stream -> Modulator -> Channel -> Equalization -> Demodulator -> Received Bit Stream

Note that I suppressed some blocks like pulse shapping, matched filter, FEC, and others. This is the most simple scheme. The transmitter modulate the signal, and sends it through a channel. In the receiver, the samples needs to be equalized in order to perform the detection and demodulation properly.

So, you can conclude that the CMA blind Equalization is somewhere between the channel and the detection.

In your code I only saw the test for the Walsh Code. I am not seeing the complete system implemented (i didn't see any modulator, for example), but I assume that now you can figure out where to include the CMA equalizer.

Let me know if you need more help.


Create an input so that poles show up as outputs? Let's say we have a rational, causal, stable LTI system with transfer function
$$H(z) = \frac{A(z)}{B(z)}$$
If $H(z)$ has $N$ poles, we can in theory have only 1 of those poles, $p_i$, show up at the output. So $y[n]$ is $p_i^n$ after some $n>A$. To do this we need to come up with an input $x[n]$ that is causal and has finite length, and $A$, so that we get that output $y[n]$.

How would we construct such a finite $x[n]$ to get this behavior?

Edit:
We can assume also that we know the most general form of $h[n]$ since $H(z)$ is causal, LTI and stable.
$$h[n] = a_0\delta[n] + a_1\delta[n-1]+a_2\delta[n-2]...a_A\delta[n-A]\\+b_1p_1^nu[n]+b_2\cdot n \cdot p_1^nu[n] + b_3 \cdot n^2 \cdot p_1^nu[n] + ...\\+c_1p_2^nu[n]+c_2 \cdot n \cdot p_2^nu[n] + c_3 \cdot n^2 \cdot p_2^nu[n] + ...\\.\\.\\.$$
So would there be a clever construction of a finite $x[n]$ to cancel all other poles except $p_i$? The first part of $h[n]$ gives us the finite first part of $y[n]$ until integer $A$.

You should always be able to achieve this after $N-1$ samples. For this I will denote the given transfer function as

$$
H(z) = \frac{A(z)}{\prod_{k=1}^N (z - p_k)}. \tag{1}
$$

Now by starting with an impulse and filter it with

$$
F(z) = \frac{\prod_{k\neq i} (z - p_k)}{z^{N-1}} = \sum_{k=0}^{N-1} \alpha_k\,z^{-k} \tag{2}
$$

and use the output of this filter as the input to the system. The combined transfer function looks like

$$
H(z)\,F(z) = \frac{A(z)}{(z - p_i)\,z^{N-1}}. \tag{3}
$$

Therefore this filtered signal should only excite the mode of the system associated with the pole $p_i$. And the $z^{N-1}$ term only adds a delay of $N-1$ samples, so it only takes a finite amount of time after which only $p_i^n$ is visible in the output.

However it can be noted that this does require you to have perfect knowledge of the poles you cancel.


Using Doppler velocimetry with digital signals such as Bluetooth I have found several sources describing how to use Doppler velocimetry/anemometry to measure the Doppler shift with analog signals, but i cannot find much about using this method for digital signals such as Bluetooth.  I was wondering if it were possible to use Doppler velocimetry/anemometry with digital signals? If so i was wondering if I could use it with the Bluetooth capabilities(2400 MHz – 2483.5 MHz) of an iPhone?

If your question is whether you could ever use the Bluetooth radio in an iPhone to perform doppler velocimetry, I'm going to go out on a limb here and say no.

Most mobile devices are incredibly  modular. This means that the BT/WiFi radio (they are often combined into a single chip these days) is designed to handle all of the RF signal processing necessary, and then hand over only the digital data received, but never the analog signals.  By the time the signal is turned into digital data, any effects of multipath, channel noise, doppler shift, and other analog phenomena has been removed.  It is highly unlikely that any mass-produced cell phone would have access to the analog signal being received by the individual radios.


Am i performing Sinesweep correctly? I am bit unsure whether i am doing i sine sweep correctly.
I am feeding my system a sine wave with f hertz within a given time interval and log its output.

And then i do the same for another frequency.

I tried this technique but seem to do it incorrect or so i think.
I tested for 2 frequencies, and fed them into  Ident in matlab and got a fit rate of 98 % which makes no sense that it could build a such a good bode plot based on 2 values..

I am doing something wrong.. should it not have more than one 2 input or is it enough?

The theory behind sweep-sine measurements of LTI systems requires a signal with constantly changing the frequency. You cannot simply playback few tones - the whole frequency range is necessary.

So that if you want to identify your system with the impulse response $h[n]$, you feed the sweep sine signal $s[n]$ into it and record the output. Obviously output will be given by the convolution:

$$y[n]=h[n]\star s[n]$$

In order to obtain the impulse response of the system, it is required to convolve the output $y[n]$ with the inverse filter $f[n]$. Generally, the inverse filter is a sweep signal, inverted in time and with scaled amplitude. So in the end:

$$h[n]=y[n]\star f[n]$$

There are two major types of sweep signals: linear and logarithmic (exponential) frequency change. The latter one is more commonly used, since its spectrum is similar to a pink noise. Personally I always use the logarithmic sweep sine:

$$s(t) = \sin\left[ \dfrac{2\pi f_1 T}{\ln\left( \dfrac{2\pi f_2}{2\pi f_1} \right)}\left(e^{\dfrac{t}{T}\ln\left(\dfrac{2 \pi f_2}{2\pi f_1} \right)} -1\right) \right]$$

where:

$f_1$,$f_2$ - initial and final frequency of a sweep

$T$ - sweep duration in seconds

There is no strict rule on what should be the length of a sweep. For example in acoustics, the rule of thumb is that it should be at least as long as predicted reverberation time.

For more literature on the topic, please refer to the following sources:

  Meng Q., et al. - Impulse Response Measurement With Sine Sweeps and
  Amplitude Modulation
  Schemes

  Farina A. - Advancements in impulse response measurements by sine
  sweeps

  Farina A. - Simultaneous Measurement of Impulse Response and Distortion with a Swept-Sine Technique


How to model Tape Saturation (Audio DSP)? I'm looking for info that would help me to build a Tape Saturation process.

I'm working with the WebAudio API. The API provides a number of DSP nodes that perform basic processes, such as a Convolver node and a WaveShaper node. The WaveShaper seems to suggest it is designed for distortion effects (http://www.w3.org/TR/webaudio/#WaveShaperNode-section). But tape saturation would seem to me more than a simple mapping of input to output values.

What are the principal processes I should be modelling? Can anyone point me in the right direction?

Saturation can be modeled to first order as a symmetrical distortion of the signal, but it also has history dependence--hence the term hysteresis in describing magnetic response to a field that is applied.  Bias (an AC field applied to the tape at the same time as the signal to be recorded) linearizes the curve and reduces hysteresis effects by eradicating the memory of the medium.  If your objective is to simulate the sound of a tape recorder, you should consider hysteresis effects, and the fact that high frequency signals are self-biasing, so that they effectively saturate at a different (lower) level than low-frequency ones.  (This was the insight behind Dolby HX processing, which really made cassette tapes sound a lot better.)


Design and implementation of causal band-pass filter for biosignals - what to consider? I currently work in a project of clinical software (C#) which deals with clinical biosignals, specifically EMG and other low-frequency signals (load-cells, goniometers, and other biomechanical instruments).
We have to be able to:

Create filters based on cuttoff frequencies (low and high) and filter order, for a given sampling frequency; The band-pass could be a cascading of one low-pass and one high-pass filters;
Apply filters in real-time, that is, adding a new input sample and immediately taking a new output sample, in a loop;

I have read a lot about filter types, but explanations are usually based on math notation, while I would need a more computational explanation.
Also, reading source code (mostly from Octave-Forge) I can only "see" a lot of numbercrunching without actually understanding what is happening.
As I understand (correct me if I am wrong), an implementation of what I need would involve:

Choose some input parameters for filter creation, that is, Low Frequency, High Frequency, Sampling Rate, and Filter Order;
Execute some magic operation on these parameters, that would give me a recurrence relation; This recurrence relation would be contained in an instance of a filter;
This filter instance would have a buffer of size N. As I start to collect data, the following loop would start to run:

Add new input sample to filter;
Apply recurrence relation;
Read new output sample;

So what I have to ask is:

Does my computational rationale about filter design and usage (as explained in the bullet-list just above) makes sense from a DSP perspective? That is, would it be a feasible plan to implement some Filter class according to such rationale?

Which concepts should I know / read about in order to properly implement such filter? For example, I understand that recurrence relation, order, sampling rate, (implicit) impulse response, are a bare minimum to define a filter, but I also read about number of taps, poles and zeros, gain, direct form, cascaded sections, and I am not sure it is necessary to deal with all the concepts to satisfy my requirements.

Is there any conceptual topic or resource (link, book) I should take a look? Preferrably one focusing on implementation (that is, without heavy math notation)?

Your bullet list makes sense and you seem to have the right idea. I think you mostly want a link that will show you how to get started, so let me suggest this.

First, start with an FIR filter. A very good FIR filter can be generated with as little as 10 lines of code, and implemented with another 10 lines of code. Basically, with FIR, if it's complicated, you are doing it wrong.

This link has source code for Windowed FIR, Parks McClellan FIR, and IIR filters, as well as code for implementing both FIR and IIR filters. Start with the Windowed FIR. They work great and are the simplest to generate.

http://www.iowahills.com/A7ExampleCodePage.html


What is the difference between the PSD and the Power Spectrum? Can any one explain the difference in the amplitude between the two and what does each one represents?

For example, when using the pwelch algorithm in MATLAB how do the two different options ('psd','power') affect the units of the output?

(I know the default is 'psd' and I can use 'power' to show something closer to the rms calculation, as stated here)

The power spectrum is a general term that describes the distribution of power contained in a signal as a function of frequency. From this perspective, we can have a power spectrum that is defined over a discrete set of frequencies (applicable for infinite length periodic signals) or we can have a power spectrum that is defined as a continuous function of frequency (applicable for infinite length aperiodic signals).

In the first case, each discrete component has units of power (W, mW, etc.). In the second case, each point in the continuous spectrum has units of power per frequency (W/Hz, mW/Hz, etc.). In the second case, one must integrate over a band of frequencies to obtain units of power. It is meaningless to talk about the amount of power at a frequency in the second case, one must talk about the amount of power contained in a spectral band (an interval of frequencies) or in several such bands. It is meaningful however to talk about the amount of power spectral density (PSD) at a particular frequency. This is an indicator of how much weight this frequency will contribute to the overall power if included in one of the spectral bands.

This allows us to compare the distribution of power at various frequencies in a signal (as opposed to comparing the power directly). This distribution is defined over a small delta of frequency (delta going to zero in the calculus sense). An alternative to PSD, would be to divide up the signal into chunks of a finite size (say 10 Hz) and compute the power for each bin (one way to obtain this would be to integrate the PSD over each bin). If another person performed a similar experiment and generated their own plot using a larger bin size (say 20 Hz), then the overall shape of both plots would be similar, but the latter plot would seem coarser and its numerical values would be greater. This would make it difficult to compare the two signals. Using power spectral density, the amplitude and bin size ambiguity is eliminated. This allows a fair comparison between
the power distribution of two different signals.

For more information, here is a link to a related Wikipedia article.


What does pre-emphasis $\alpha=1$ will do to a signal? I am trying to understand pre-emphasis equation.

$$ y[n]=x[n]-\alpha x[n-1]$$

Where $\alpha = 0.95$.
My question is what does it mean if $\alpha = 1$?

The pre-emphasis filter is a simple high-pass filter with a monotonically increasing magnitude of its frequency response $H(\omega)$. Its value at DC is $H(0)=1-\alpha$, and its value at Nyquist (i.e., the maximum frequency) is $H(\pi)=1+\alpha$. Consequently, for $\alpha=1$ you get a zero at DC, which means that the filter will totally suppress a constant input signal (which is obvious from looking at the difference equation).

Note that a pre-emphasis with $\alpha=1$ cannot be compensated for because the DC component is lost and can't be recovered. In signal processing terms: there is no stable inverse filter.


Difference: Ergodic channel capacity and spectral efficiency Can someone explain what is the difference between the capacity and ergodic capacity? Ergodic capacity and spectral efficiency?

Loosely speaking:

Capacity is the supremum of data rate that one can send data with arbitrarily small error probability over a given channel;
ergodic capacity is also the supremum of rate with arbitrarily small error probability, but for fading channels, in assuming the fading process is ergodic, as the term "ergodic" suggests;
spectral efficiency is data rate divided by occupied bandwidth.

Therefore,

Ergodic capacity is not capacity. Capacity is not ergodic capacity. They are channel-specific and we do not necessarily know how to achieve them.
Spectral efficiency is system-specific, i.e. we have already had a system, as said MarcusMuller in his comment. We ignore error probability in using this. Nonetheless, a spectral efficiency is usually implicitly linked with an "acceptable" error probability.

In short, they are very different performance indicators with very different assumptions. In academy, you are free to choose any indicator that supports your research.


What methods can be used to compute homography I am trying to understand panorama stitching. I've come to a conclusion that homography is one of the fundamental parts in the process. Now I am trying to make sense of all the information from literature across the web concerning about solving for homography.

I've read about Direct Linear Transform, Least Squares method and RANSAC. However, from what I've read I'd say DLT and RANSAC are the methods to solve homography while Least Squares is just used as additional step in DLT when there are more than 4 correspondence points. Is this right? I feel like I have things mixed up.

I've never actually done Panorama Stitching, but I was talking to people who were, and I was involved in a project where we had to estimate the homography for fine localization for visual navigation, so here it goes.

Basic approach

Rough outline of the approach using RANSAC (you can read more of the approach in this Master thesis, page 12,13 - unfortunately, in Croatian):

estimate pairs of corresponding points for two images $A$ and $B$ (lot of them) and get rid of the obviously bad correspondences

(detailed summary of a feature matching approach in my Master Thesis, pages 13-20 - In English)
pick any four points for the initial estimate of homography $H$

(most details on the numerical methods used to estimate the homography can be probably be found in my Bachelor Thesis (in Croatian), but there should also be some details in this Master Thesis (Cro) and this one (Eng))
project all points from the image $A$ onto the image $B$ using the first estimation of the homography $H$
compare the projected points from $A$ to the actual points present in $B$ (using Euclidian distance)
discard all the point pairs where the distance between projected point from $A$ and corresponding point from $B$ is greater than a threshold $t$ -- those points are outliers
(this can be repeated several times, possibly with a more constrictive threshold $t$ every time, to get rid of more outliers)
the points that are left are our population points
using only the population points, do the final estimation of homography $H_{final}$ just from those points

Some additional info

Coincidentally, if you have the initial rough matching results, the rest of this is implemented in OpenCV, it's a function called findHomography and there is a good StackOverflow post on that. It includes RANSAC (and from what I can read from that post, Least Squares method should have a similar goal as RANSAC, but RANSAC gives better results).

If you need better estimates, you can find references for methods other than RANSAC in this Master thesis (Eng) - pure RANSAC approach is compared to something called MLESAC, and RANSAC coupled with Least Median of Squares method.

Suggestions on an efficient implementation of Homography estimate from points (presuming that you do not want to use OpenCV or some other available implementation) can be found in my Bachelor Thesis (Cro).

To estimate the initial (not-very-precise) pairs of correspondences, you should use matching methods. Those methods are based on:

detecting feature points

(there's a lot about that in this DSP question, and especially in this answer you will find a link to the current widely used framework for comparing feature detectors)
using appropriate feature descriptors

(nice info in this DSP question)
doing matching between those points and possibly some kind of filtering to further better the matches obtained

Matching a point from an image $A$, $x_{A}$ with a descriptor $d(x_{A})$ is usually done just by finding the point in the image $B$, $x_{B}$ with the most similar descriptor $d(x_{B})$ according to Euclidian distance.

The most simple (and surprisingly efficient) filtering method is to compare the distance from $d(x_{A})$ to the closest descriptor $d(x_{B_1})$ and second closest descriptor $d(x_{B_2})$. If $\frac{D_{eucl}(d(x_{A}), d(x_{B_1}))}{D_{eucl}(d(x_{A}),d(x_{B_2})}$ is greater than some coefficient $k \in [0,1]$ (typically $k = 0.6$), the match is rejected. First introduced by D.G. Lowe: Distinctive image features from scale-invariant keypoints.

Additional filtering can be done according to spatial arrangement of the corresponding points. One such interesting method was described in J. Sivic, A. Zisserman: Video Google: A Text Retrieval Approach to Object Matching in Videos.


Protect an IIR filter from being reverse-engineered I created a somewhat unique IIR filter and I want to protect the filter from being reverse-engineered

I think you all know it is quite easy to get all the different weights of an IIR by using impulse signals.

The filter is programatically encapsulated in a bigger program, so I can add some tricks to protect the filter coefficients.

But what would be the best method to protect my IIR filter ?

Short answer:

You can't. If an attacker can insert a signal that covers the whole bandwidth (e.g. a white signal, or at least one that has no spectral zeros) into the system (and he can do that over an arbitrarily long time, or add up observations), they will get an output, and can through the magic of correlation get the impulse response.


Procedural wave generation, frequency shift issue sorry for the long post but I want to try and get all the details in here in hopes it will help.

I've been working on building a procedural audio system for use in XNA projects. First this is general structure of the procedural system.

WaveFragments: generate the actual wave forms for a specified duration, based on frequency using a phase accumulator.

Waves: complete wave forms, generated by joining wave fragments together.

DynamicWaves: use a wave fragment to generate small buffers of waveform data, which are pushed into a DynamicSoundEffectInstance to be played.

My problem is (I think) with my wave generation. At first it seems to work OK, but then after a certain time the sounds frequency alters slightly, and then back again later.
for example a tone of C3 (130.81f) will increase frequency slightly (but noticeably) at around 20 seconds and then back again later on at 77 seconds. (I've not yet had the patience to see if it shifts again later, I've gone as far as 3 minutes without it shifting up again.

I'm going to assume its a issue with the way I'm generating the wave data with my phase accumulator.

Here is my generation of the PhaseDelta based on provided frequency.

 _PhaseDelta = (float)((Math.PI * 2 * _Frequency) / (float)samplesPerSecond);

and here is the calculation I use for the generation of the wave.

for (int i = 0; i < numSamples; i++)
{
     _PhaseAccumulator = MathHelper.Lerp(_PhaseAccumulator, _PhaseAccumulator+_PhaseDelta, _PhaseSmoothness);
     data[i] = (float)Math.Sin(_PhaseAccumulator) * _Amplitude;
}

Thanks in advance for any input that anyone might have. It may be the most glaringly obvious fix, but at this point after hours of dull single tones my brain is two melted to be getting anywhere with fixing this now.

Now I've tried it without lerping the accumulator (obviously produces the clicks at the start) but still produces the shift. I've also tried lerping the phase delta gradually to its new value, again to no avail.
 (EDIT: i have also tried throwing lerp out of the window and just adding the phase delta each sample, again no dice).

So, my question is, any ideas on why this weird pitch shift is happening? In practicality it isn't a huge issue, the system was only ever intended for short notes to be stacked together for old video game style music, or a tone that would likely change frequency often. But I'd feel more comfortable knowing that the wave generation is accurate and wasn't accumulating some error somewhere sending it way of track.

I suspect that the accumulator value is growing very large, which causes the float's resolution to degrade.  To avoid this you need to occasionally reduce the accumulator to manageable levels.  One way to do that is to subtract by 2*pi whenever you exceed 2*pi.


Removing white noise by taking the mean of many samples I'm learning about signal processing, and I am attempting to remove white noise from an output of a blackbox system.

As you can see from the above image (left is the input sin(t), right is the output) there is quite a bit of noise.
I am attempting to verify that it is white noise. To do, I have run the system 100 times with the same sin(t) input and have taking the mean of the outputs, which resulted in the following plot:

If it is white noise, I believe I should expect that the mean of many outputs will approximate the noise-free output for larger and larger sample sizes.
However, as you can see, there is still quite a bit of noise and I'm not too sure why. At first I thought the data between samples wasn't lining up when taking the mean, but I don't believe that should actually matter.

The variance of the white noise will go down by $N$ through averaging (or standard deviation will go down by $1/\sqrt{N}$, where the standard deviation as a magnitude quantity will be visually consistent with the "spread" of noise on the graphic the OP is looking at. Thus to make it visually reduce by a factor of 10, the OP would need to average 100 times (if it is indeed white noise).
In practice with actual signals there is usually an upper limit on how much we can average, until non-stationary noise contributions (such as phase noise, 1/f noise etc) limit further improvement. If further averaging does not reduce this further, then such noise processes may be limiting further improvement.
Also be very careful in a simulation of white noise that the random number generator isn't actually creating the exact same sequence each time (in which there will be no change through averaging)- given the OP's plots I suspect that is what may be occurring here. This would be clear by inspecting each sample in successive runs. The other confirmation of a white noise process is to take compute the autocorrelation function  on the waveform, which for the white noise component will be an impulse at $\tau=0$ (or experimental approximation of one, limited by the actual time duration of the signal, so a large spike proportional to duration of the signal, and very small everywhere else).
One approach to remove white noise from a black box system with an unknown signal confined to a specific band of frequencies is to pass the output of the signal through a bandpass filter which performs a form of weighted averaging over previous samples of the one signal (rather than what the OP has done which is essentially averaging a system with N outputs, which by definition of ergodicity would be equivalent if the system is ergodic.)


What kind of filter would have a sawtooth frequency response? I want a complex FIR or IIR filter whose response is approximately a sawtooth in the frequency domain. How should this be formally described, and what formula or algorithm can I use to calculate taps for one?
Here is a sketch of the magnitude response I want (not analysis of an actual filter). I think the phase response doesn't matter.

A windowed design is easy. Starting with an ideal zero-phase frequency response:

$$\frac{\omega}{2\pi} + \frac{1}{2}$$

The inverse Fourier transform of that gives the ideal impulse response:

$$\int^\pi_{-\pi}\left(\frac{\omega}{2\pi}+\frac{1}{2}\right)e^{i\omega k}d\omega = \frac{\sin(\pi k)}{k} + i\left(\frac{\sin(\pi k)}{\pi k^2} - \frac{\cos(\pi k)}{k}\right).$$

Similar to the $\text{sinc}$ function, the real part $\frac{\sin(\pi k)}{k}$ has an illegal division by zero at $k=0$ where the limiting value $\pi$ should be used instead. For other integer $k$, the real part is zero. The imaginary part has similar problems at $k=0$, where its limit is 0. Drawing from Rick's answer, at integer $k$, $\sin(\pi k)$ is zero and $\cos(\pi k)$ is +1 or -1. This simplifies the ideal impulse response at integer $k$ to:

$$\begin{cases}\pi&\text{if }k=0,
\\i\frac{-(-1)^k}{k}&\text{otherwise.}\end{cases}$$

If the filter is implemented as a sum of a purely real and a purely imaginary filter, the purely real filter requires no work other than to multiply the input by a constant. The imaginary part of the impulse response is antisymmetrical, which may enable savings.

To continue with the windowed design, choose a window function and its length (equals the number of filter taps) and multiply the ideal impulse response by it.
That gives the impulse response of the filter. The window function should be centered at $k = 0$.


System that has derivative of input is non causal Consider a system $y(t) = \dot{x}(t)$ where $y$ is the output and $x$ is the input. Given an initial condition $x_0$ and two inputs $x_1$ and $x_2$ such that $$x_1(t)=x_2(t) ,  0 \le t < t_0$$ the outputs $y_1,y_2$ corresponding to $x_1,x_2$ respectively should satisfy $$y_1(t)=y_2(t) ,  0 \le t < t_0$$ since $$x_1(t)=x_2(t) ,  0 \le t < t_0 \implies \dot{x}_1(t)=\dot{x}_2(t) ,  0 < t < t_0$$
How is this system non causal then?

I am not sure why what you wrote is relevant for the causality.
Anyways, the derivative operator is depends on the future signal.
Actually this almost causal, in the sense that you need an infinitesimal environment around t to define the derivative.

Just for example, we can take 2 cases:
$$ x_1(t) = -t , x_2=|t| \Rightarrow \dot x_1(t)=-1, \dot x_2(t) =undefined$$
$$\Rightarrow y_1(t)\neq y_2(t)$$
despite $\forall t\leq 0 : x_1(t)=x_2(t)$


Revert anayltic signal(length(N/2)) to real signal (length(N)) As a beginner in signal processing, i'll try to explain my problem most thoroughly.
I'll firstly explain the idea of the transformation the real valued data has undergone.
This includes the ideas and reasons why it does work and the way to revert the procedure. Then i'll explain my difficulties in reverting it.

First of all a property of the FFT of a real signal:

Given an $N$ data point real valued signal $x(t)$, it has $\frac{N}{2} + 1$ unique Fourrier coefficients $X(f_i)$ with  $i \in {0,..,N-1}$. The other $\frac{N}{2}-1$ coefficients can be retrieved by inverting the spectrum of the coefficients $X(f_1)$ to $X(f_{\frac{N}{2}})$.
So in some sense all the information of my signal is stored in the first $\frac{N}{2}+1$ of it Fourrier coefficients. Keep that in mind, because this property i'm going to exploit to retrieve the original signal from it's analytic representation.

Now to the theory of the procedure the data has undergone:

Let $X_a (f_i)$ be the analytic representation of the Fourrier coefficients $X(f_i)$ of $x(t)$.
Then $X_a(f_i)=0$,  if $f_i < 0$, and for the Fourrier Coefficients of the original signal $x(t)$ we know that $X(f_i)= \frac{X_a(F_i)}{2}$with $i \in 1,...\frac{N}{2}$. So $\frac{N}{2}-1$ coefficients carry no additional information about the signal. This allows a downsampling to $\frac{N}{2}+1$ frequencies, by cutting those $X_a(f_i)$ off, with $f_i <0$. And no information is lost, as the original signal was real. To make sure you can follow me http://en.wikipedia.org/wiki/Analytic_signal explains it in detail. Let's call the resulting complex signal of length $\frac{N}{2} +1$ $x_c(t)$.

The way to retrieve the original real signal would be:

If i'd want to retrieve my original signal all i'd have to do is invert the Fourrier Cofficient vector $[X_a(f_i)]$, $i \in 1,N-1$ and append them to my frequency spectrum.
Resulting in the spectrum: $[X_a(f_0),...X_a(f_{\frac{N}{2}+1}), X_a(f_{\frac{N}{2}}),X_a(f_{\frac{N}{2}-1}), X_a(f_{\frac{N}{2}-2}),...X_a(f_1)]=FFT(x(t))$ of my original signal $x(t)$.

Now to my situation:

I do get something similar to the results of above modification of signals.
The major difference is that i don't get all $\frac{N}{2}+1$ coefficients of the complex signal $x_c(t)$. But what i do get is the complex signal from $0$ to $\frac{N}{2}-1$. So what i'm lacking is the coefficient $X_a(f_{\frac{N}{2}})= \sum_{k=0}^{N} e^{i\pi k}x_a(k)$ of the Fourrier transform of my Fourrier pair $(X_a(f_i),x_a(t))$. This is a real valued Fouerrier coefficient of the signal $x(t)$ at the nyquist frequency.

So what way i would be able to get the best approximation of my original signal?

Thanks in advance!!

PS: Hope this format is way better than before.

First answer deleted since I completely misunderstood the question.

The problem here is that apparently you are missing one frequency point. From the question it's not clear what is missing. You say you have the points from 0 to N/2. That's actually N/2+1 (if 0 and N/2 are included) and all you need. If that's not the case, there are gew possibilities

Go to the person who provided the data and ask them to fix it
X(0) and X(N/2) are actually real valued. Sometimes people pack the X(N/2) value into the imaginary part of X(0) to save space and create a symmetrical representation. Check your X(0) and figure out whether it's imaginary part is actually 0 or not
There is a way to represent a N-point real sequence with exactly N/2 complex coefficients that involves a  half sample rotation in the frequency domain. That's pretty esoteric and unlikely in this case (but works quite well).
If the high frequency energy in your original signal is sufficiently low, you can set the Nyquist coefficients simply to zero. According to the sampling theorem, that should be zero anyway. In practice this will add some noise to the reconstructed signal but if it's small enough that may be tolerable.


Wireless Body Area Networks with Minimum Energy Consumption For adaptive compressive sensing(cs),the sensing matrix is related to the input signal.
For example, in rakeness-based(cs), the sensing matrix is obtained by solving an optimization problem which maximizes the correlation between the input signal and the sensing matrix,meanwhile minimizing the coherence between the columns of the sensing matrix.
I know that the operations at the sensing node should be minimized, to minimize the energy consumption.

Where is this optimization problem performed ? At the sensing node ?? If so, how can this rakeness- based approach implemented in a wireless body area network without consuming much energy??

There are two optimisation problems here:

To get the sensing matrix $\mathcal{A}$
To recover the signal by $y = \mathcal{A} x$

Out of these two problems, the first one is the one that is solved off-line to obtain a sensing matrix $\mathcal{A}$ that appears to be more "suitable" to a particular family of signals.

The second optimisation problem is the standard compressive sensing recovery and this is the bit that might take more CPU effort to be solved.

For an indication about the kind of effort required, you might want to have a look at this page, the code it provides and specifically the Optimization/l1eq_pd.m file.

This is the entry point to a generic CS "solver" that then calls other optimisation algorithms contained in Optimization/.

If you take a look in there you can see that you can expect matrix multiplications, matrix inversion, point to point multiplications of sequences and other less intensive operations.

All of these also depend on the amount of data captured, so you might have to do some experimentation with realistic signals from the domain you are dealing with to get a more accurate idea of how much computational effort is required and therefore how much power to solve these problems.

It certainly looks more like a problem for a Raspberry Pi Zero rather than an Arduino type of CPU, if you are considering a mobile / wearable application.

Hope this helps.


How to shift a signal on the $y$-axis to match another signal in MATLAB I have two signals recorded in the time-domain then converted to the frequency domain. These signals are $z$-axis accelerometer measurements recorded from an Arduino 101.

Now visually they are fairly similar (with the peaks and falls sort of lining up), however the difference is in their amplitude, e.g.:

I will be comparing the power of the two signals, so I need them to be in-line with each other, ie. shift either one of the signals to fall around the same area.

I'm not too sure what the terminology for this is, but any help would seriously be appreciated.

Matlab's polyfit function can be used to do linear regression. This can be used to find a linear relationship (offset and gain) that maps one dataset to the other.


DCT and Hard Thresholding If I have an Image and i find the DCT and then apply hard thresholding on the coefficients and then IDCT then I have attenuated the noise. Can someone please explain in detail or point me to the answer as to why this works. I understand why a filter blocking high frequency will work in denoising(because we assume that noise consists of high frequency components) but why does magnitude thresholding work?

Low-pass filtering is based on the assumption that "natural" images have more energy in the low-frequency coefficients than in the high frequency coefficients ; while noise will affect equally all coefficients. Thus, removing the high frequency coefficients will comparatively eliminate more noise than signal. The problem is that there are "legit" high-frequency coefficients in images, such as edges. Denoising by low-pass filtering will attenuate those coefficients and cause edge blurring.

The property used in magnitude thresholding is sparsity - "natural" images are very likely to have only a small set of high non-zero coefficients in the frequency domain. Adding uniform (and independent from pixel to pixel) noise is equivalent to adding a small random value to all the frequency coefficients. The result is that all the coefficients that were 0 in the original image now have a small value ; while the coefficients that were high in the original image are comparatively untouched. Thresholding the small magnitude coefficients will cancel the contribution on noise on those coefficients - though it will not recover the effect of noise on the high coefficients.

If you replace edges by transients/attacks, the same thing applies to sounds by the way.


Scaling the input vs scaling the impulse response for an LTI system Two different cases:

We pass $x(t)$ to an LTI system with impulse response $h(2t)$ and get the output $y(t)$.

We pass $x(2t)$ to an LTI system with impulse response $h(t)$ and get the output $z(t)$.

I'm looking for possible restrictions (e.g. on the input / impulse response) such that the $y(t)$ and $z(t)$ are equal, if possible.
So far I have this:
$$
y(t) = x(t) * h(2t) \implies Y(j\omega) = X(j\omega) \frac{1}{2} H(j\omega/2)
$$
$$
z(t) = x(2t) * h(t) \implies Z(j\omega) = \frac{1}{2} X(j\omega/2)H(j\omega)
$$
If we take, for example, $X(j\omega)$ from the first equation and plug into the other one, we get
$$
Z(j\omega) = \frac{Y(j\omega/2)H(j\omega)}{H(j\omega/4)}
$$
I'm stuck at this point to comment on the relation between $y(t)$ and $z(t)$ and assert some restrictions to make them equal possibly.

Aside from

Input that is only a constant value or
a System that does nothing to the signal, i.e. $H = \text{const.}$,

these things are not the same.
Sketch of a proof:

Scaling input shifts energy to frequency components that hadn't energy before,
Scaling the system means changes the Transfer Function, but keeps the system LTI, which means it can't add energy to any components that had none before.


Apply filter before or after Downsampling? Context: I'm trying to do formant estimation. My first attempt didn't work out all that well so I'm trying to follow closely Praat since their formant estimation is remarkably accurate.
Question: I thought it was obvious to apply the filter (pre-emphasis) before Downsampling, however Praat doesn't. Is there a reason for this? Are filters usually applied before or after downsampling?

I haven't looked at the reference yet , but just to answer the question, yes filtering is performed prior to down sampling to avoid aliasing back in the out of band signals. However if your signal is already appropriately bandlimited you could forego the filter entirely. There's various reasons you might want to do other, unrelated, filtering after the downsampling operation. In fact, doing such filtering after the downsampler is preferable due to the lower sample rate and thus lower computational complexity.


FFT returns large low-frequency power - probably because signal is interpreted as containing cycles that do not fit. How is this called? So I have a large number of signals like the one in the first picture below and I would like to extract and compare the frequencies within them.
I applied a Fourier transform which resulted in results like those shown in the second picture. The large power of very low frequencies suggests to me that the signals include something that is interpreted as being of such low frequency that one cycle does not fit into the time window (which is 1.2 seconds).
I hence applied a High pass filter with different cut-offs and the pattern is gone after using a HPF with 2Hz (last picture).

I have two questions about this:
1) Does my reasoning make sense?
2) How is such an effect / pattern called?

I need to include this in a paper I am currently drafting and hence am especially unsure about question 2). My supervisor's first reaction was that these are edge-effects, but looking into them I feel like that is not really what is at work here.

Thank you!!

Your reasoning sounds correct. It is impossible to see from your example what frequency that low frequency is.

If the frequency of the phenomenon is 0 Hz, then it is simply a direct current (DC) bias that affects the first bin of the discrete Fourier transform (DFT). If none of your signal of interest is in that bin then you can simply ignore it in your analyses or set it to zero. If your signal never goes to negative values, then you have DC bias (although there could be other low frequencies).


How to derive equation for second order interpolation of soft knee cutoff in a compressor? I am reading Audio effects: theory, implementation and application by Joshua D. Reiss and Andrew McPherson.
They explain that, in the simplest case with just a hard cutoff threshold, you have a piecewise linear equation
$$y_G = \begin{cases}
x_G & x_G <T\\
\frac{x_G -T}{R} + T & x_G \geq T
\end{cases}
$$
with a corresponding diagram like

However, in many cases, you may not want this hard cutoff. You may want it to be softer like this

To do this, Reiss proposes the following formula

I found a PDF of a presentation he did where he called this "second order interpolation", but I have no idea what that means other than I think it's second order because of the quadratic terms in this version.
How does Reiss derive the middle formula in this piecewise equation?
Reference:

https://www.eecs.qmul.ac.uk/~josh/documents/2012/GiannoulisMassbergReiss-dynamicrangecompression-JAES2012.pdf
http://c4dm.eecs.qmul.ac.uk/audioengineering/compressors/documents/Reiss-Tutorialondynamicrangecompression.pdf
Differences between soft knee and hard knee in Dynamic Range Compression (DRC)?

If you look at the formula and plot it, you get what you see in the other link on dsp.ee. But if you add the inflexion points and plot the consituent functions, it should become more clear:

The function would have been a linear one, with two slopes, but between the points defined by $T-W/2$ and $T+W/2$, there is a curve uniting the two segments. This means that the conditions become:

a simple, quadratic formula (to avoid too much computation)
the function's value at the first point should be equal to the point, itself, and
the derivatives at those points should be the actual slopes

With these, there are many ways to interpolate, so I can't say exactly how they did, but these must have been their starting points.
For example, considering two slopes of $3/4$ and $1/4$ at $[-1,1]$, the quadratic function $f(x)=ax^2+bx+c$ that would unite those would have to have:
$$\left\{
\begin{aligned}
f'(-1)&=2a(-1)+b&=\dfrac34 \\
f'(1)&=2a(1)+b&=\dfrac14
\end{aligned}\tag{1}
\right.$$
The derivative makes $c$ disappear, but its value at $-1$ is $3/4$, so the solving the equations above results in the function:
$$f(x)=\dfrac34+\dfrac{x}{2}-\dfrac{x^2}{8}\tag{2}$$
and this is how it looks like:

The slopes at $[-1,1]$ are a match, which means the two slopes are now connected by a quadratic formula. The whole graph would be, from left to right: the red slope until $-1$, then the quadratic blue until $1$ and then the green slope; just like in the 1st picture.

I'll add that you can go the other way, just as well: from the derivative, up.
You know that the signal has a slope in the beginning, and one at the end which means that their derivatives give you fixed values (in the example above they were $3/4$ and $1/4$). Uniting the points will be a slope of the form $ax+b$, as seen in (1), starting from $-1$ with the value $3/4$, and ending up at $1$ with $1/4$. Since it's a line, it will pass through the middle point, $m$, and it will have a slope $s$, both given by:
$$\begin{align}
m=\dfrac{\dfrac14+\dfrac34}{2}&=\dfrac12 \\
s=\dfrac{\dfrac14-\dfrac34}{2}&=-\dfrac14 \\
\Rightarrow\quad g(x)=m+sx&=\dfrac12-\dfrac14x\tag{3} \\
\end{align}$$

Integrate this to get the quadratic, to which you add the constant term $3/4$ -- the value of the function at $-1$:
$$\int_x{g(x)\mathrm{d}x}=\color{blue}{\dfrac34}+\dfrac{x}{2}-\dfrac{x^2}{8}\tag{4}$$
(2) and (4) are the same. q.e.d.


LSM level set method vs thresholding I am trying to understand the LSM algorithm applied to grayscale image segmentation.

There are essentially 2 things that are blocking me:

1) From my point of view - the level sets - i.e. ''moving'' the 3D surface (that is, conceptually visualizing the grayscale image (say 8 bits/pixel) as a surface whose height is given at each point (pixel $x,y$) by the intensity value of that pixel, through a fixed plane, or the other way around, is exactly the same as performing a thresholding from 0 to 255 (i.e. each "level set" would be all the pixels at threshold value 0,1,2 ...255)?

2) Let $f(x,y)$ be the image function (associating an intensity between 0 and 255 to each point $(x,y)$. The LSM
is defined as $$\Gamma = \{(x,y) | \phi(x,y)=0\}$$

However, it does not make sense to me, (from my point of view which is probably wrong, hence my question), that $$\phi(x,y)=0$$ as as the surface has its own height at each point and it would make more sense to me to consider that the plane intersecting the surface is moving from bottom to up, i.e. I can ''accept'' better the formulation given for the level sets:

$$\Gamma = \{(x,y) | \phi(x,y)=c\}$$

and for the example I use, $$c \in [0,255] $$

Question 2 is more like a ''conceptual'' problem, but question 1 is more fundamental, since I cannot seem to see any difference with a simple thresholding.

Rephrased:

Question 1) Can we consider (loosely speaking, as indeed as the li is a set of points not a matrix) that each level curve $$l_i$$ of a level set of an image I is (in programmatic terms, Matlab syntax): li= (I==i) ? (where i is the height, in grayscale 8 bits, i is a number between 1 and 255, and as Matlab uses matrices everywhere, I is a $n \times m$ matrix and the operation == checks for every pixel in I if it is equal to i and the result is a binary image with 0s or 1s, i.e. this is called thresholding, hence my question is LSM the same as thresholding with operator == (and the single '=' is the assignment operation of course))

Question 2) Is it essentially (conceptually) the same to imagine moving the shape through a fixed plane, or moving a plane through a fixed 3D shape ?

What I was (in an incorrect way) trying to say is that we can visualize the level curves $$l_i$$ as the white circles on the thresholded image (although not perfect since it is discretized an maybe rounding errors/noise etc but suppose they are continuous) and ideed i do understand that the result of thresholding is a matrix so strictly speaking not equal to a lvl curve. But "loosely" speaking, visualizing the lvl curves do you agree?

Question 1) Can we consider that each level curve
  $$li$$
  of a level set of an image I is [...] a binary image with 0s or 1s, i.e. this is called thresholding, hence my question is LSM the same as thresholding with operator == (and the single '=' is the assignment operation of course)).

No it is not.

$\Gamma$ is a curve. The result of thresholding is a matrix.

$\phi$ is positive inside the boundary, negative outside of the boundary and $0$ right on the boundary.

A function that is a good candidate for $\phi$ would be the distance function. An alternative formulation to that, if you really wanted to represent the boundary of a closed region would be to obtain a slice via thresholding and then apply a mask such as:

$$\begin{matrix} 1 & 1 & 1 \\ 1 & -8 & 1 \\ 1 & 1 & 1 \end{matrix}$$

This would give you the pixels that are right on the boundary of the thresholded region. You still have to extract them of course and if we are talking about a curve you would also have to connect the points (otherwise a set of pixels in 2D space is not necessarily a curve) and while connecting them you would also have to select a winding (clockwise or counter-clockwise). In addition, the thresholding might return regions that are not connected in which case, you would have to decide if they are part of the same object and (optionally) use morphological operators to make sure that they are closed curves and not paths.

  Question 2) Is it essentially (conceptually) the same to imagine moving the shape through a fixed plane, or moving a plane through a fixed 3D shape ?

While I do understand where this question comes from, I have to say that in this particular case it is neither, because $\phi$ is not the image.

The meaning of the level set as "all points for which $\phi$ is zero" is not the same as sliding a threshold over $\phi$ which is what you are trying to denote with the second formula.

Yes, you do need to slide the threshold so that you segment regions based on the way they "fill up the basins" of an image. But for each threshold you need to calculate a separate $\phi$.

The reason for this is because two successive $\phi$s (e.g. a given slice and a slice below) will inform your boundary tracking about which way to move towards, or in the level set, which pixels to try to mark as belonging to the same region. This is because, once you have your $\phi$ you then put a derivative on it and the derivative is the "instrument" that informs you which way the boundary is going (expanding, contracting, remaining the same).

Hope this helps.


Passband Attenuation I just started out in signal processing I am trying to make a low pass filter with a 0 - 10000Hz Passband and has a maximum attenuation of 3db. Is the passband attenuation the same as the passband ripple? If so why is it that even though I am using a Bessel filter I still have a passband ripple of -3db?

Passband attenuation refers to the maximum signal loss within the passband of the filter which should ideally be zero; i.e., pass band shall not attenuate the signal it passes; it's a unit gain region.
The term ripple is about of alternating gain through the passband for certain filter types yet it's still used even if the passband is not alternating such as a Butterworth filter. In that case the frequency at which the gain reduces to -3dB is the cutoff frequency for a lowpass filter. Or similarly high and low cutoff frequencies for a bandpass filter.


reconstructing time-stretched signals for filtering purposes I would like to know if it's possible to reconstruct the original time domain signal from it's time-stretched version? Is there any algorithm out there that can do this? Python, Matlab, etc?

I want to isolate certain sound events in the signal, however, the events are really close together, that it's hard to filter. A solution that I had in mind was to stretch the signal, filter this time-stretched signal and then reconstruct back by back stretching using inverse stretching factor to get the original length.

The algorithm will of course have to keep the correct pitch at all times. Is this possible or not?

Keeping your “correct pitch” after a time stretch can become a more difficult problem than your original problem of filtering short events.  Both require a time-frequency resolution trade-off.


number of multiplications per second required by the system Consider a system implementing a rational sampling rate change by 5/7: for this, we cascade upsampler by 5, a lowpass filter with cutoff frequency π/7 and a downsampler by 7. The lowpass filter is a 4rd-order Butterworth filter with transfer function

$$H(z) = \frac{b0 + b1z^{-1} + b2z^{-2} + b3z^{-3} + + b4z^{-4}}{1 - a1z^{-1} - a2z^{-2} - a3z^{-3} - a4z^{-4}}$$
Assume that the input works at a rate of 1000 samples per second. What is the number of multiplications per second required by the system? Assume that multiplications by zero do not count.

Since there are 9 filter coefficients and 5000 samples left (after upsampling, filtering, downsampling), 9*5000 = 45000 is my MAC calculation, but that seems to be wrong approach, any other idea to solve this question???

One potential option would be to split the lowpass filter into the FIR part and the pureley recursive IIR part.
$$H(z) = H_1(z) \cdot H_2(z) = \frac{1}{1 - a1z^{-1} - a2z^{-2} - a3z^{-3} - a4z^{-4}} \cdot \left [ b0 + b1z^{-1} + b2z^{-2} + b3z^{-3} + + b4z^{-4} \right ]$$
You still have to apply the recursive part to each sample in the upsampled domain, but you only need run the FIR when you need an actual output sample. So this could be
$$N = 5000*4 + (5000*5)/7$$
That's not a great idea in terms of stability and noise performance, but would work if you have enough precision and head room.
UPDATE:
There are further ways you could potentially reduce multiplactions. The FIR section of a butterworth can be split into two second order sections which are simply $[1 2 1]$ so you could potentiall write this as
$$H(z) = H_3(z) \cdot H_4(z) \cdot H_4(z) = \frac{b_0}{1 - a1z^{-1} - a2z^{-2} - a3z^{-3} - a4z^{-4}} \cdot \left [ 1 + z^{-1} + z^{-1} + z^{-2}  \right ] \cdot \left [ 1 + z^{-1} + z^{-1} + z^{-2}  \right ]$$
So in this case it's down to 25000 multiplies per second, although it's not a particularly useful implementation. At the end there are multiple different ways of implementing them and which one is better depends on the properties of the platform.


Design the coefficients of Wide-band FIR filters It will be a great help if someone could share a content/tutorial about how to design the coefficients of Wide-band FIR Filters.

As it stands your question appears to be quite broad, so the answer depends largely on your application and level of knowledge.

If you have no prior experience, you might want to consider picking up a book about digital signal processing, e.g.

Digital Signal Processing by John G. Proakis and Dimitry Manolakis

which may be available in your local library.

For digital filters specifically, there are great resources available online such as this free book:

INTRODUCTION TO DIGITAL FILTERS WITH AUDIO APPLICATIONS by JULIUS O. SMITH III

There is also the Cookbook for audio EQ biquad filters by RBJ which you should also check out, found for example here (although i cannot find the original link right now).


Frequency Hopping Spread Spectrum Synchronization Methods Frequency-hopping spread spectrum waveforms change the carrier frequency at some interval.

I can see if the frequency is constant for the duration of an entire packet, the pilot symbols could be used for synchronization but what happens in the case of "burst-by-burst carrier frequency hopping"? As mentioned in this paper, it is possible to change the carrier frequency within a symbol duration.

In order for the receiver to perform coherent demodulation (like for PSK), how would phase synchronization occur in the case where the carrier frequency changes at the symbol rate (ie each symbol is transmitted at a different carrier frequency)?

This can be accomplished by changing the carrier frequency using a Numerically Controlled Oscillator (NCO) which maintains an accurate and continuous phase versus time trajectory via the phase accumulator. This is markedly different than changing the frequency with a classical PLL where we would typically break and reacquire lock to change frequency resulting in a random starting phase position at each carrier. Since the desired phase location is repeatably known at each symbol starting position, we can measure phase error to correct for carrier offset.

This is further demonstrated with the following graphics:

A typical non-hopped coherent receiver with a carrier offset would have a phase error versus time such as that shown in the plot below (where the specific modulation is not shown).

Since the question is not how to do carrier tracking in general, I assume for purposes of this post that approaches to recover and track the carrier in this circumstance are well understood.

Now consider a frequency hopped system, where we have used an NCO to change the carrier frequency. Without a carrier offset (the signal as observed at the transmitter), the phase versus time for three successive symbols, with each at a new carrier frequency, would appear as follows (again with the specific modulation used not shown). Frequency is the derivative of phase, so a constant frequency offset result in a linear slope of phase, as shown.

With the same carrier offset as we first presented, the phase trajectory would actually be as follows in the receiver:

Even though we are changing frequency, because the phase is continuous and predictable (no random phase inserted at the start of each symbol as in the PLL switched case), we are still able to determine the background phase slope (error) that is used to drive the carrier tracking loop.

For more details on NCO see: Numerically Controlled Oscillator (NCO) for phasor implementation?

For more details on Carrier Tracking see:
High modulation index PSK - carrier recovery


First Order Hold discrete-time approximation to first order continuous-time linear system Consider the following first order linear system described by:
$\tau \frac{d}{dt}y(t)=-y(t) + x(t)$.
I have seen a discrete time approximation to this system using a "First Order Hold" approximation of the form:
$y[k+1] = e^{-T/\tau} y[k] + [-e^{-T/\tau} + \tau/T(1-e^{-T/\tau})]x[k] + [1 - \tau/T(1-e^{-T/\tau})]x[k+1]$ where $T$ is the sampling period.
Can anyone show a derivation of this result?

The solution comes from assuming that
$$
x(t) = x[k] + \frac{x[k+1]-x[k]}{T}(t-kT), \ \ t\in[kT,(k+1)T)
$$
where $x[k]:=x(kT)$, this is, a first order hold or a linear interpolation of $x[k]$ and $x[k+1]$ in between $kT$ and $(k+1)T$. Now, we simply compute the EXPLICIT solution of the differential equation (at instants $t=kT$) by integrating it from $t=kT$ towards $t=(k+1)T$. This should be enough to derive the result but let me get into more detail. First, multiply the ODE by $\frac{1}{\tau}e^{t/\tau}$ to obtain:
$$
\frac{dy}{dt}e^{t/\tau} + \frac{1}{\tau}e^{t/\tau}y(t) = \frac{1}{\tau}e^{t/\tau}x(t)
$$
and notice that $\frac{d}{dt}\left(e^{t/\tau}y(t)\right) = \frac{dy}{dt}e^{t/\tau} + \frac{1}{\tau}e^{t/\tau}y(t)$ and thus:
$$
\frac{d}{dt}\left(e^{t/\tau}y(t)\right) = \frac{1}{\tau}e^{t/\tau}x(t)
$$
Now integrate both sides:
$$
\int_{kT}^{(k+1)T}\frac{d}{dt}\left(e^{t/\tau}y(t)\right)dt = \int_{kT}^{(k+1)T}\frac{1}{\tau}e^{t/\tau}x(t)dt
$$
the left hand side results (by fundamental theorem of calculus) into:
$$
\int_{kT}^{(k+1)T}\frac{d}{dt}\left(e^{t/\tau}y(t)\right)dt = \left.e^{t/\tau}y(t)\right|_{kT}^{(k+1)T}  =e^{(k+1)T/\tau}y[k+1] - e^{kT/\tau}y[k]
$$
Thus,
$$
e^{(k+1)T/\tau}y[k+1] = e^{kT/\tau}y[k] + \int_{kT}^{(k+1)T}\frac{1}{\tau}e^{t/\tau}x(t)dt
$$
Now, simply divide by $e^{(k+1)T/\tau}$ to obtain
$$
y[k+1] = e^{-T/\tau}y[k] + e^{-(k+1)T/\tau}\int_{kT}^{(k+1)T}\frac{1}{\tau}e^{t/\tau}x(t)dt
$$
As you can see we are almost done. We just have to compute:
$$
\begin{aligned}
&e^{-(k+1)T/\tau}\int_{kT}^{(k+1)T}\frac{1}{\tau}e^{t/\tau}x(t)dt =\\ &e^{-(k+1)T/\tau}\int_{kT}^{(k+1)T}\frac{1}{\tau}e^{t/\tau}\left(x[k] + \frac{x[k+1]-x[k]}{T}(t-kT)\right)dt
\end{aligned}
$$
which we can do explicitly since $x[k]$ and $x[k+1]$ are just constants. After some algebra you can factor terms with $x[k]$ and terms with $x[k+1]$ and you will obtain precisely what you wanted. If you want to check, I plugged this into wolfram alpha here ("x1" is $x[k]$ and "x2" is $x[k+1]$).
This kind of discretizations are nice since, if $x(t)$ does behave as a piece-wise linear continuous function, the discretization is NOT an approximation, but an equivalent (and exact) representation of the ODE, this is: since it was obtained by explicit computation of the solution of the ODE, no information is lost about the output samples $y[k]$. (in contrast to approximating $dy/dt \approx (y[k+1]-y[k])/T$ or some other approximation as the bilinear one)
NEW REMARK: I just want to point out that this same procedure can be performed for other high order systems, or even for other types of holders. Recall that any LTI system with scalar input $x$ and scalar output $y$ can be represented by its transfer function, or equivalently by a state space representation
$$
\frac{dz}{dt} = Az + Bx,\ \ \  y = Cz
$$
with $z(t)\in\mathbb{R}^n$ with $n$ the order of the system, $A$ a square matrix, $B$ a column vector and $C$ a row vector. Recall that the poles of the transfer function are the eigenvalues of $A$ and zeros of the transfer function depend only on the coefficients on contained in $C$. Thus, we can solve for the system on $y[k]$ explicitly using a similar procedure as before:
$$
z[k+1] = e^{AT}z[k] + e^{A(k+1)T}\int_{kT}^{(k+1)T}e^{-A\tau}Bx(\tau)d\tau, \ \ y[k]=Cz[k]
$$
where $e^{A\tau}$ is the exponential matrix. Finally, plugging in the particular form of $x(t)$ results in the final form of the discretization (in discrete time state space representation). This can be used to derive (exact or "consistent") discretizations for ZOH, FOH, high order holders or more exotic holders.


Z transform - Inverse System function - Why number of poles and zeros myst be equal? I know that if a system is causal then the system function H(z) must have :
a) a ROC that spans from the exterior of the most distant pole and
b) the number of zeros must not be greater than the number of poles

I have found this exercise:

and the solution:

Why the number of poles and zeros must be equal?

Regardless of causality and stability, if you count poles and zeros at the origin and at infinity, the total number of poles always equals the total number of zeros. I'll show a few examples to make this obvious.

First, take a polynomial

$$P(z)=(z-z_0)(z-z_1)\ldots (z-z_M)\tag{1}$$

It has $M+1$ zeros and no (finite) poles. However, it must have $M+1$ poles at $z=\infty$ (because of the term $z^{M+1}$).

Next consider an "all-pole" function

$$A(z)=\frac{1}{(z-z_0)(z-z_1)\ldots (z-z_N)}\tag{2}$$

with $N+1$ poles and no (finite) zeros. Clearly, due to the $N+1$ poles, the function $A(z)$ in $(2)$ has $N+1$ zeros at infinity, because of the term $z^{N+1}$ in the denominator.

In general, if you have

$$H(z)=\frac{(z-z_{0,0})(z-z_{0,1})\ldots (z-z_{0,M})}{(z-z_{\infty,0})(z-z_{\infty,1})\ldots (z-z_{\infty,N})}\tag{3}$$

and if $M>N$, you get $M-N$ poles at infinity. For $M<N$, you have $N-M$ zeros at infinity. Consequently, in all cases the number of poles equals the number of zeros if poles and zeros at infinity are included.


Estimate sampling frequency given an array of values Given an array x of size 2000, is there any way to calculate (or estimate) the sampling frequency using Matlab? The values are the only information I have about the signal. I know the signal has been sampled with Fs = 360 Hz but I want to know if there is any way to extract this information from these plots.

This is the signal plot, it is an ECG:

And this is the FFT:

Given an array x of size 2000, is there any way to calculate (or estimate) the sampling frequency [...]? The values are the only information I have about the signal.

No.

  I know the signal has been sampled with Fs = 360 Hz but I want to know if there is any way to extract this information from these plots.

No.

The sampling frequency is an intrinsic characteristic of a discrete sampled signal. Any operation on the sampled signal uses Fs either implicitly or explicitly. The only way to infer the sampling frequency is to inject or detect a pilot tone on a properly sampled signal.

For example, the Discrete Fourier Transform does not make explicit reference of the sampling frequency. The DFT simply decomposes the signal into the sum of those sinusoids that can be supported by the "data frame" of length $N$ that it is presented with. But, when you then want to go from "wave number" or "bin", back to the physical frequency it represents, then you have to know the Fs.

If you don't use the "right" Fs, then you might run into the condition of aliasing where you basically cannot trust the numbers you are looking at because they are an uncontrolled overlap of the signal spectrum and itself. This is why, the anti-aliasing filter is a key component in the analog pathway of the signal BEFORE the analog to digital conversion. The anti-aliasing filter will make sure that no frequency very near or above $\frac{Fs}{2}$ will make it through the system.

If we were to remove this filter, the digital processing part's outputs would become unreliable.

Having said this and as is pointed out by @Fat32, if you are absolutely sure that the sampled signal adheres to this restriction and it was sampled without aliasing then you can infer the Fs. BUT, this cannot be done with a random Electrocardiogram. Or rather, it really depends on the sort of accuracy you are after.

A standard electrocardiogram waveform's harmonics can extend to 100Hz.

The question now is, what is the actual bandwidth of THIS specific ECG that you have in your hands? That one, may indeed have been obtained with a standard montage but it contains noise and baseline wandering. Also, depending on what the patient was doing, the signal would be different and its bandwidth would be different too.

Therefore, which DFT bin corresponds to that 100Hz is uncertain.

If you are interested in a ball park figure, you can look at a point in the signal's DFT where its power drops below an average background noise and take that point as your 100Hz nominal value and work backwards to what Fs was used. But that would certainly not be accurate, especially when "strange" Fs are used (like 128.5 Hz).

So, in the general case, it is impossible to "second guess" the sampling frequency from a discrete signal's values.

The only exception to this is to ensure that no aliasing will occur (i.e. sample high enough) and use a pilot tone. That is, a tone or other signal whose spectrum you know exactly or can control fully.

If there was a single pilot tone in the signal then you could say, I know that the pilot tone is at 1kHz, it shows up on the DFT output as a high output on bin k, then I can work backwards to estimate the correspondence of bins to physical frequencies.

Could there be such a pilot tone in your signals?

Yes, it is likely that your signals contain power line interference. This will show up as a consistent set of spikes in the DFT because it is not only that you would be looking at a sinusoid of 50Hz but you would get a few harmonics along too.

The frequency of the power line is not exactly 50Hz (in Europe at least) but it has to adhere to strict standards, so perhaps for the relatively short duration of the ECG it could be used as a pilot tone at an assumed value of 50Hz.

Of course, it is not necessary that the ECG contains power-line noise. It may have been suppressed by a filter or other shielding at the point of acquisition. In that case, you are back into "it cannot be done" response.

Hope this helps.


Butterworth high pass filter My signal's base line is around 1. When I apply Butterworth high pass filter with cut off frequency Fc = 2Hz signal's base line drops to 0. My data are from the accelerometer. Is this behavior normal ?

Yes.  A high-pass filter removes low frequencies, and a constant baseline is the lowest frequency possible ("DC" = 0 Hz), so it will be completely removed and the baseline will become 0.


Solving LASSO ($ {L}_{1} $ Regularized Least Squares) with Gradient Descent To the best of my knowledge, state of the art methods for optimizing the LASSO objective function include the LARS algorithm and proximal gradient methods.

I was wondering however, if the LASSO objective function
$$ ||y-Ax||_2^2 + \lambda ||x||_1$$ can be optimized using (vanilla) gradient descent?

I know that the 1-norm is not differentiable at zero, but numerically one is never exactly at zero. Therefore the gradient exists at each iteration point $x^{(k)}$ and I do not really need to resort to subgradients.

Unlike for the L2-penalized least squares, the gradient descent step would of course need a decreasing step size and would not yield exact zeros.

Can this method be called gradient descent or would it be better to call it subgradient descent or something else?

Many thanks!

It can easily solved by the Gradient Descent Framework with one adjustment in order to take care of the $ {L}_{1} $ norm term.

Since the $ {L}_{1} $ norm isn't smooth you need to use the concept of Sub Gradient / Sub Derivative.
When you integrate Sub Gradient instead of Gradient into the Gradient Descent Method it becomes the Sub Gradient Method.

In the case of the $ {L}_{1} $ norm the Sub Gradient is given by the function $ \operatorname{sign} \left( \cdot \right) $.
In this context it works on vector in Element Wise manner.

So the Pseudo Code (MATLAB) would look something like:

for ii = 1:numIterations
    vG = (mA.' * (mA * vX - vY)) + (paramLambda * sign(vX)); %<! Sub Gradient
    vX = vX - (stepSize * vG); %<! Sub Gradient Update
end

Yet, As mentioned by @David, this is not a modern way to solve this problem.
Iteration wise it seems that Coordinate Descent is the fastest method for solving the LASSO ($ {L}_{1} $ Regularized Least Square) problem.

You may look at a project I created (${L}_{1}$ Regularized Least Squares - Solvers Analysis) which compares many method for that optimization problem:

Implementation is included so you can see how it works.


Infinite extent of spectrum, but also in time in Oppenheim's Discrete Time Signal Processing? In Oppenheim's Discrete Time Signal Processing there's on p. 323 no limited band in both time and frequency - wouldn't that violate the Heisenberg Principle?

Not at all. The Uncertainty Principle says that a function cannot be both limited in time and limited in frequency. More specifically, the product of the signal's widths in time and in frequency (i.e., its time extension $\Delta_t$ and its bandwidth $\Delta_f$) is bounded from below:

$$\Delta_t\cdot\Delta_f\ge C\tag{1}$$

where the constant $C$ depends on the definition of bandwidth and time extension.

Note that $(1)$ is a lower bound, not an upper bound, so both widths can be infinite without contradicting $(1)$.

If a function is sharply localized in time then, by the Uncertainty Principle, it cannot be sharply localized in frequency, and vice versa. However, if - as in your example - a function is NOT localized in one of the two domains then this does not mean that it must be localized in the other domain. In may very well be non-localized in both domains, as is the case in the given example for non-integer $\alpha$.

Also take a look at this question and its answers for more details on the the Uncertainty Principle.


Calculate aliasing of $x_a(t) = \cos{(2\pi300t)} + \cos(2\pi600t)$ when sampled with $F_s = 1000$ I'm asked to sample the signal
$$x_a(t) = \cos{(2\pi300t)} + \cos(2\pi600t)$$
with sampling frequency $F_s = 1000$ and plot the magnitude spectrum for the resulting sampled signal.

My thinking is that the frequency of $300$ does not change, so it just results in the normalized frequency $\frac{3}{10} = 0.3$.

Since $600$ is above the nyquist frequency, aliasing arises. So we get the frequencies $\frac{6}{10} = 0.6$ and $\frac{6}{10} - 1 = - 0.4.$

So in total I would plot peaks at the frequencies $\pm0.3, \pm0.4, \pm0.6$, but the answer is supposedly:

Why? Where is my thinking wrong? And what is the method to always get the right aliasing peaks?

Peter said it all in his comment, but I'll try to explain it here in a bit more detail.

Your thinking is correct, but incomplete. Of course, the component at $f=300$ (unit depends on normalization of $t$) will not cause aliasing, so there will be no additional component below Nyquist, but what you forget is that sampling makes the spectrum periodic with period $f_s$ (sampling frequency). So the $f=\pm 300$ component results in components at $nf_s\pm 300$, $n\in\mathbb{Z}$. For $n=1$ you get the components at $f=1000\pm 300$, one of which (the one at $f=700$) is the one you see in the figure at normalized frequency $0.7$.


Basic Fourier Series Not Understanding how to interpret n values I have a question regarding $n$ in the Fourier series. If a question states "find the Fourier series for (any function) and find the values of an $b_n$ etc", say I'm finding it from the trig way when $n<=4$. What $n$ values do I use, is it 1,2,3,4 and if so why does it start at 1 and not include 0. Also why doesn't the $n$ values keep going past negative so e.g 4 3 2 1 0 -1 .... What makes it start at 1 for this value of n.

Is it since in the formula n=1,2,3.....?

However if it stated for $-3<=N<=3$, would the n values be -3 -2 -1 0 1 2 3? if so why? Also would the zero be in it or not? My $n$ values were less than 0 would you then have to use the complex way?

I'm assuming the "n" you are referring to is the coefficient index, also referred to as the harmonic number; i.e. the n=3 term refers to the "third harmonic".

The n=0 term is the average value of the function, and is sometimes called the "DC component" because it does not depend on time.

Depending on the version of the Fourier series you are using, the terms in the series could go from negative infinity to positive infinity, or they could just go from n = 0 to positive infinity. The trigonometric version of the Fourier series (with sine and cosine) starts at n=0 and goes to infinity, whereas the complex exponential version (i.e. $c_n e^{j\omega_0 n t}$) has terms for negative and positive values of n.

The two versions are equivalent because the negative n terms cancel out the imaginary parts of the positive n terms in the complex exponential series.


Series Sum Operator Properties and Identities As I was reading my book, I came across the following expression:
\begin{equation}
\sum _{n=-N}^N\sum _{m=-N}^N \phi \left(n-m\right)
\end{equation}
Then the book said, if we let $k = n-m$, then:
\begin{equation}
\sum _{k=-2N}^{2N}\left(1-\frac{\left|k\right|}{2N+1}\right)\phi \left(k\right)
\end{equation}
Can somebody please tell me if this is a property of summation operator and how they got there. Also please refer me to a reference that carries such properties of summation. Thank you in advance.

That equation is almost correct. Tim gave the correct solution in a comment, which is just a scaled version of the given formula. I'd like to show how you can quite easily see this.
First, note that
$$\sum_{m=-N}^N\phi(n-m)=(r\star\phi)(n)\tag{1}$$
where $\star$ denotes convolution, and $r(n)$ is a sequence of ones in the range $n\in[-N,N]$. Second, summing a sequence is also a convolution with the sequence $r(n)$ evaluated at $n=0$. So we can write
$$\begin{align}\sum_{n=-N}^N\sum_{m=-N}^N\phi(n-m)&=(r\star r\star\phi)(n)\Big|_{n=0}\\&=(s\star\phi)(n)\Big|_{n=0}\tag{2}\end{align}$$
where $s(n)$ is the convolution of two rectangular sequences, which is a triangular sequence:
$$s(n)=(2N+1)\left[1-\frac{|n|}{2N+1}\right],\qquad n\in[-2N,2N]\tag{3}$$
Combining $(2)$ and $(3)$ gives the desired result, which equals the formula in your question, scaled by $2N+1$.


How do I find the Hilbert transform (qudrature) of a discrete signal by convolution with a Hilbert kernel? Given a discrete real signal $ f_{re}(t) $ the (complex) analytic signal is given by:
$ x(t) = f_{re}(t) + j*f_{im}(t) $.
I want to calculate $f_{im}(t)$: the quadrature by convolving with a Hilbert kernel of size 2*n+1:
$f_{im}(t) = H*f_{re}(t)$
$H=[H_{-n},H_{-n+1},..,H_0,H_{n-1},H_n]$
According to "Complex seismic trace analysis",M. T. Taner*, F. Koehler*, and R. E. Sheriff;
$ H_i = 2*sin^2(\pi i/2)/(\pi i), H_0 = 0$
for $ n = \infty $.
I want to use a finite length kernel, e.g. n=16.
Do I need to taper the kernel from the equation above with a window?
How do I normalize the odd (antisymmetric) result kernel?
Am I right in assuming that calculating the quadrature by convolving with H introduces a smoothing?
If so; should I also smooth the real signal, e.g. by a gaussian kernel of size 2n+1?
EDIT:
In Matlab there is a function called "firpm" that creates a "Parks-McClellan optimal equiripple FIR filter".
It seems this is generally considered the "optimal" FIR approximation to an IIR?
Does this mean that using the simple coefficient equation from "Complex seismic trace analysis" is too simplistic?
EDIT2:
regarding the normalization: I tested in Matlab and found that convolving with a kernel calculated from the simple coefficient equation gave good results in the case of the input signal being a sine. However I had to divide by a factor m to get correct normalization:
H5 = [-0.0000   -0.6366  0    0.6366    0.0000]: m = 0.1271
H7 = [-0.2122 -0.0000 -0.6366  0  0.6366 0.0000 0.2122]: m=0.2525

What is the general formulae for m?
EDIT 3:
The ideal kernel would have constant magnitude response. However it seems that since my kernel will be of finite size; it will act as an bandpass filter, and therefore I need to look out for bandpass ripple.
The solution seems to be a magnitude response that rolls of gently at the low and high frequencies. The "simple coefficient equation" is quite good in this respect; it does not introduce a lot of ripple, instead the shorter the kernel the earlier it starts rolling off. As an example here is the 19 point kernel:

In comparison the 19 point kernel from Matlabs firpm does not roll of at high frequencies (!?) but otherwise have roughly the same amount of ripple as my simple kernel:

My conclusion so far thus seem to be: no the "simple coefficient equation" gives a good approximation to the optimal FIR for the Hilbert convolution operator. The "Parks-McClellan optimal equiripple FIR filter" might be slightly more "optimal" but the difference is neglible, at least for reasonably large kernels (say 19 points and upwards). Anything I am forgetting? I will therefore move on with this approach. Only remaining problem: how do I normalize?

The short answer: with a finite convolution you can't calculate the Hilbert transform exactly.

The long answer: You will have to make approximations and decide what types and amounts of error you can live with. This will depend heavily on the application, the signals itself and what you are trying to do with the results. This paper may give good guidance on math and potential trade offs http://www.andrewduncan.ws/air/index.html

Per se, the Hilbert Transform is not a smoothing operation. For example the HT of a cosine is a sine, which isn't any smoother.


Online materials on Support Vector Machine I am doing my final year project on Image Processing.I want to know about the Support Vector Machine.Please suggest some links for the former.

Thanks and regards in Advance

Video Lectures : http://nptel.iitm.ac.in/courses/117108048/

Notes : http://www.nptel.iitm.ac.in/courses/106108057/ (Read from Lesson 25. These notes helped me very much)


Which codec and format are Audible downloads? Files downloaded with the Audible Download Manager are stored with the .aax file extension, but the information I can find about this format seem to vary.
Regarding this source, the files are stored in the M4B container (MPEG-4 Part 14), also a quick look into the file header of the .aax files reveals aax M4B mp42isom.
Furthermore, regarding this source, the audio is encoded with the ACELP codec.
If I would decrypt the .aax files DRM with my Audible authentication code without converting the audio, which format and codec would the audio be in?
(Also I'm not sure if this belongs here or in "Sound Design".)

As this post over at the apple stackexchange says, .aax files from Audible are indeed encrypted m4b files, which means that if I would decode an .aax file without converting it, I would get a .m4b file.
This is possible with the following command:
ffmpeg.exe -activation_bytes XXXXXX -i input.aax -c copy output.m4b

The -c copy ensures that the data will not be converted.


Who is Barbara (test image) In recent years I've started doing work in image analysis and have used the Barbara and Lena images. Early on I became aware of the interesting history surrounding the Lena image and that got me wondering: Who is Barbara, where was the image shot, and what was the original intent (presumably it wasn't created for image research)?

Many people cite Allen Gersho's lab or Marco Schmidt's (now defunct) test image database. However, this is a bit unsatisfying to me.

The appeal of this image is obviously in the numerous lines, which test the aliasing properties of resizing, denoising, and super-resolution algorithms.

It seems Allen Gersho is the source, according to the Acknowledgement section of Embedded image coding using zerotrees of wavelet coefficients

If you want more specific information you could ask him or one of his students from the early 90s, when this image first seems to have appeared.


When sinusoidal input starts at n=0, why are transient response associated with z-transform poles of digital filter? In http://www.eng.ucy.ac.cy/cpitris/courses/ECE623/presentations/DSP-LECT-10-11-12.pdf, it says that when sinusoidal input $X(z)$ starts at n=0 (with n<0 having zero input) and the input passes through digitl filter $H(z)$, resulting in output z-transform $Y(z) = H(z)X(z)$, transient response is associated with poles of the filter while steady-state response is associated with the poles of the input. Is this true? If so, how does one prove this?

Yes, that's true if the $\mathcal{Z}$-transform of the input signal has poles on the unit circle (i.e., the input is a constant or a sinusoid, starting at $n=0$), and if the filter is stable, i.e. all its poles are inside the unit circle. If you split up $Y(z)$ by a partial fraction expansion, you get a contribution for each pole. The contributions from the poles of $H(z)$ all decay to zero because we've assumed that the filter is stable, i.e. all its pole magnitudes are smaller than $1$, and consequently, these terms decay exponentially. On the other hand, the contributions from the poles of $X(z)$ do not decay because the corresponding poles lie on the unit circle. This means that after a while, all contributions associated with the poles of $H(z)$ will be very small (practically zero; transients), whereas the contributions from the poles of $X(z)$ remain at a constant amplitude (steady-state response).


Why doesn't law of large numbers apply to this stationary time-series? There's a paragraph in Wikipedia that states the following:

  Let Y be any scalar random variable, and define a time-series $\{X_t\}$, by
  $$X_{t}=Y\qquad {\text{ for all }}t$$
  Then $\{X_t\}$ is a stationary time series, for which realisations consist of a series of constant values, with a different constant value for each realisation. A law of large numbers does not apply on this case, as the limiting value of an average from a single realisation takes the random value determined by $Y$, rather than taking the expected value of $Y$.

I cannot understand the reason for the law of large numbers not being applicable. Any explanation would help.

I believe that you are thinking that each value of $X_t$ is determined by a different realisation of $Y$, which in this example is not true.

Suppose that $Y$ is the value that comes out from a dice throw. Thus, $Y$ can achieve any integer value between $1$ and $6$. Throw the dice once. Suppose that you get $2$. Then, according to the definition in the example, your time-series would be:

$$X_t=2 \qquad {\text{for all }}t$$

What I think the quote from Wikipedia is saying is that the mean value of this stochastic process doesn't coincide with the mean value of the random variable $Y$, as

$$\mathbb{E}[X_t]=2 \neq 3.5 = \mathbb{E}[Y]$$

If you constructed the time-series letting each value of $X_t$ be determined by a new throw of the dice, then the means would match.


Do convolution encoders exist with two separate (independent) inputs? Do convolution encoders exist with two separate (independent) inputs?

I got the question in my mind when I saw the question attached here.

Do convolution encoders exist with two separate (independent) inputs?

Yes, such convolutional encoders do exist, and the diagram that is included in your question is an example of such an encoder. The code that is generated is called a (3,2) convolutional code.

Now what is your real question?


How to calculate sine wave SNR from oscilloscope I'm trying to calculate the SNR of a sine wave produced with a 12 bit DAC. I'm new when it comes to calculating/measuring noise in signals. My idea was to capture the signal using an oscilloscope (which unfortunately does not have an FFT function), so I could run an FFT analysis in MATLAB to separate the signal in its fundamental and harmonics to find the SNR.
This is the signal that I acquired, it's a 600 mVpp sine wave at 1kHz

I have 625000 samples, and according to the metadata from the oscilloscope, it has a sample interval of 1.6e-8s. My first approach was to use MATLAB's snr function directly using the voltage samples
Fs = 625000;
r = snr(measurements, Fs);

Using that code I obtain this:

I was expecting the fundamental to be close to 1kHz, but it's sitting at 0.011 kHz with -9.699 dB.
Is this result expected?
I also tried plotting the FFT using the following code (which again gives me this peak close to 0):
y = fft(measurements);
n = length(measurements);
fshift = (-n/2:n/2-1)*(625000/n);
yshift = fftshift(y);
plot(fshift,abs(yshift))

For all the reasons above, an oscilloscope is inadequate to this task.  Here's a couple of alternatives:
Method 1: find someone with a really good spectrum analyzer, and borrow time on it.
Method 2: build a really good notch filter (for a one-off you'll want to use analog components, and good ones).  Characterize it (I didn't say this would be easy).  Then run your 1kHz signal through the notch and look at the remainder in the O-scope.  Once you get rid of the 1kHz tone, whatever remains is the stuff you don't want -- characterize that, add the 1kHz tone back in on paper, and you've characterized the net effect of your DAC and your notch filter.
Note that method 2 requires some significant analog chops.  If your go-to solution involves white proto-board, junkbox capacitors and LM741 amps, you probably need to up your game.


Finding the memory polynomial model coefficients I want to build a memory polynomial model given that I know the corresponding output signals for known input signals.
What is the easiest way to find the coefficients of it?
And how do we take the memory and order into consideration while we build it.

This is what I know. if A and B are input and output signals of length n, then A_inverse x B will give me n by n  matrix.
But how can I introduce memory and order?, if i add any extra columns to A_inverse x B, it will change the dimensions to n by m where m > n. and i can no more operate it on A.
 'x' means matrix multiplication.

Memory polynomial system identification is a linear regression problem.  Suppose you sent a $N$-element sequence $x(n)$ through a system and obtained an $N$-element output sequence $y(n)$.  If the system were described by a memory polynomial, then your model would be...

$y(n) = \sum_{q=0}^{Q}[x(n-q)a(q)+x(n-q)|x(n-q)|^2 b(q)] + \epsilon(n)$

Here $a(q)$ and $b(q)$ are coefficient of the models and $\epsilon(n)$ is an error term.  This is a linear model in $a$ and $b$, and can be expressed in matrix form as...

$\begin{bmatrix}y(Q)\\y(Q+1)\\...\\y(N-1)\end{bmatrix} = \begin{bmatrix}x(Q)&...&x(0)&x(Q)|x(Q)|^2&...&x(0)|x(0)|^2\\x(Q+1)&...&x(1)&x(Q+1)|x(Q+1)|^2&\ldots&x(1)|x(1)|^2\\\ldots&\ldots&\ldots&\ldots&\ldots&\ldots\\x(N-1)&\ldots&x(N-1-Q)&x(N-1)|x(N-1)|^2&\ldots&x(N-1-Q)|x(N-1-Q)|^2\\\end{bmatrix}\begin{bmatrix}a(0)\\...\\a(Q)\\b(0)\\...\\b(Q)\end{bmatrix}+\begin{bmatrix}\epsilon(0)\\\epsilon(1)\\...\\\epsilon(N-1)\end{bmatrix}$

You can solve for $a$ and $b$ using any of the linear least squares functions in MATLAB, Python, etc.


Estimating the amplitude of a particular frequency I am relatively new to DSP and have been reading a lot on the internet. I have a couple of questions.

I have a signal in the form of a function
$$
f(x) = A_0 + A_1 \cos(\omega_1 x) + A_2 \cos(\omega_2 x)+...+ A_n \cos(\omega_n x).
$$
I have $f(x)$ and I know the minimum and maximum frequencies. $f(x)$ can have really high frequencies.

a) I need to find the amplitude $A_p$ of a particular frequency $\omega_p$ in $f(x)$. One way to do this is, I could get sample points from $f(x)$ sampling at greater than Nyquist rate and do an FFT and find the amplitude. But if $f(x)$ has really high frequencies, the sampling would take most of the time. If i am interested only in estimating the amplitude $A_p$, is there another way, a faster way, of doing this. In the best case, I would like to estimate $A_p$ to within 0.5 of its actual value.

b)If I want $A_0$, the DC component, I can compute the average of $f(x)$ but with high frequencies in $f(x)$ that would mean high sampling rates. Any other faster way of doing this?

c)The FFT takes $O(n\log(n))$ if $n$ is the size of the input and returns the magnitudes for $n$ different frequencies. If we are only interested in one particular frequency we should be able to alter the source code so that it would take only $O(\log(n))$, right?

Thanks for answering!

The usual way to estimate the amplitude of a particular frequency is to use the Goertzel algorithm.  There is a good write-up by Rick Lyons here.

Even though Rick's writeup is about single tone detection, it can be applied when multiple tones are present, too.


How to upscale an image without affecting pixel aspect ratio? For example, if we upscale an image from 100 x 100 to 400 x 300 then we see the pixel stretch in the image which I want to eliminate.

The easiest way I could think is to crop the image to maintain the pixel aspect ratio i.e. upscaling to 300 x 300 or 400 x 400

Edit (Clarification): I want upscale an image without causing pixels to stretch, upsampling algorithms like nearest neighbor, bilinear interpolation, etc cause pixel stretching, especially when there is a change in aspect ratio e.g. from 1:1 to 16:9

Popular cheap methods are:

Letterboxing

Cropping

A combined method: cropping, overlapped by the letterboxed image

The combined method works surprisingly well in videos, and people are quite accustomed to viewing such videos.

Photoshop and others also have smart content-aware scaling/resizing algorithms that re-organize or generate image content to the new image size. I don't have PhotoShop to demonstrate it with, but here is a video from Adobe.


Identification of a transfer function Can you obtain a transfer function in frequency domain with a bode plot form given by this function:
$$ S=k_1 \sqrt{1+\frac{k_2}{f}} $$
I did not managed it because actually in the lowest part of frequencies decreases with 10dB/dec. So is not a simple pole (20dB/dec) and zero function.
Maybe it is possible to obtain such function by adding many low pass filters.
Any idea?

The objective is to produce noise, presumably Gaussian distributed, that has both white and 1/f components. This is evident from the desired power transfer function, i.e., the square of the OP's "S" expression, which would be used to convert input white noise to output white plus 1/f noise. The  OP does not state whether the PSD should be unilateral or bilateral, so it will be assumed, without loss of generality, to be unilateral. Hence the unilateral noise PSD is

$$\mathrm{S^2 = k_1^2(1+k_2/f)} \tag{1}$$

where $\mathrm{k_1^2}$ is the unilateral noise PSD of white noise only, with units of $\mathrm{V^2/Hz}$ and $\mathrm{k_2}$ is the noise corner frequency, with units of $\mathrm{Hz}$. In a digital simulation,

$$\mathrm{k_1^2} = 2\sigma^2 /f_s = \sigma^2 /f_{nyquist} = 2\sigma^2 \times \Delta t    \tag{2}$$

where $\sigma$ is the standard deviation of the Gaussian white noise (units: V), $\mathrm {f_s}$ is the constant sampling frequency (units: Hz), $\mathrm {f_{nyquist}}$ is the Nyquist frequency (units: Hz), and $\mathrm {\Delta t}$ is the constant sample time spacing (units: s).

To generate the noise, first select parameters in accordance with equation 2. For example, if it is desired to have $\mathrm{k_1^2} = 5.00 \times 10^{-8}$ $\mathrm{V^2/Hz}$ and $\mathrm {f_s} = 1$ $\mathrm{kHz}$, then $\sigma = 0.005$ $\mathrm{V}$. This is the standard deviation entered into the Gaussian white noise generator in the figure below:

The second step is to specify $\mathrm{k_2}$. Using the noise generator described previously here, it is not feasible to directly enter a noise corner frequency into the non-white noise generator code. However, it is necessary to enter the desired low limit and high limit for the non-white noise, as shown in the dialog box in the figure, so the high limit can be adjusted empirically to obtain the desired noise corner frequency. Note that $\gamma = 1$ generates $\mathrm{1/f}$ noise and there is an optional gain adjustment as well.

The next figure shows a temporal trace of the noise produced with $\mathrm{k_1^2} = 5.00 \times 10^{-8}$ $\mathrm{V^2/Hz}$, $\mathrm {f_s} = 1$ $\mathrm{kHz}$, $\sigma = 0.005$ $\mathrm{V}$ and the non-white noise parameters shown in the first figure:

Then the third figure shows the average of 1000 32k FFT PSDs:

The noise corner frequency is where $\mathrm{f = k_2}$. From the plot, this is about 1.6 Hz or so.

If this noise was the dominant one in an experiment, and if the signal was baseband, then the measurement would be afflicted with $\mathrm{1/f}$ noise and all the annoying consequences that entails. But if it was feasible to modulate production of the signal, before the signal encountered the noise, then the signal could be deliberately located in the white noise 'flats' and a lock-in amplifier could perform the customary demodulation.